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ABSTRACT 

In two previous papers l|Price fc Monaghanll2004alrJ) (papers 1,11) we have described 
an algorithm for solving the equations of Magnetohydrodynamics (MHD) using the 
Smoothed Particle Hydrodynamics (SPH) method. The algorithm uses dissipative 
terms in order to capture shocks and has been tested on a wide range of one dimen- 
sional problems in both adiabatic and isothermal MHD. In this paper we investigate 
multidimensional aspects of the algorithm, refining many of the aspects considered in 
papers I and II and paying particular attention to the code's ability to maintain the 
V • B = constraint associated with the magnetic field. In pa rticular we implemen t a 
hyperbolic divergence cleaning method recently proposed bv lDedner et alJ l)2002|) in 
combination with the consistent formulation of the MHD equations in the presence of 
non-zero magnetic divergence derived in papers I and II. Various projection methods 
for maintaining the divergence- free condition are also examined. Finally the algorithm 
is tested against a wide range of multidimensional problems used to test recent grid- 
based MHD codes. A particular finding of these tests is that in SPMHD the magnitude 
of the divergence error is dependent on the number of neighbours used to calculate 
a particle's properties and only weakly dependent on the total number of particles. 
Whilst many improvements could still be made to the algorithm, our results suggest 
that the method is ripe for application to problems of current theoretical interest, such 
as that of star formation. 

Key words: (magnetohydrodynamics) MHD - magnetic fields - methods: numerical 
- star formation 



1 INTRODUCTION 

Magnetic fields play an important, in some cases crucial, role 
in many areas of astrophysics. Despite the relative simplic- 
ity and well-studied nature of the equations which describe 
them, their effects are complicated and analytic studies are 
difficult and limited in scope. It is for this reason that a large 
theoretical effort over the past decade or so has been devoted 
to developing accurate numerical algorithms for Magnetohy- 
drodynamics (MHD) in an astrophysical context. There are, 
however, severe technical challenges to be overcome in the 
numerical solution of the MHD equations. 

Smoothed Parti cle Hydrodynamics (SPH, for a review 
SCO lMonaghan|[l99l is a fully Lagrangian particle method 
which solves the equations of fluid dynamics on a system 
of moving interpolation points which follow the fluid mo- 
tion. SPH is an extremely versatile and robust numerical 
method and as a result has found widespread use in Astro- 
physics. There have, however, been difficulties with previ- 
ous attempts to simulate magnetic fields within SPH, most 



prominently due to a numerical instability found to occur 
when an exactly momentum conserving form of the SPMHD 
equations was used. 

In two previous papers JPrice fc Monaehanl Eooiallbl 
hereafter papers I and II), we have described an algorithm 
for SPMHD in detail. The discrete equations are formulated 
from a variational principle (paper II) which ensures consis- 
tency with physical principles (such as conservation of mo- 
mentum and energy) and a consistent treatment of magnetic 
divergence terms, the effects of which we will investigate in 
this paper. Artificial dissipation terms appropriate for shock- 
type problems were formulated in paper I. These terms are 
carefully formulated to give a positive definite contribution 
to the entropy. The algorithm has been tested on a wide 
range of standard one dimensional problems used to test re- 
cent grid-based MH D codes and a l so on the one dimensional 
'Toy Stars' of lMonaehan fc Pried <20u4 . The algorithm has 
been shown to give robust and accurate results on these 
problems. 
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In more than one spatial dimension errors associated 
with the non-zero divergence of the magnetic field need to 
be taken into account in any numerical MHD scheme. There 
are two distinct issues to be addressed. The first is the treat- 
ment of terms proportional to V • B in the MHD equations 
(in particular in the formulation of the induction equation 
and the magnetic force). The second is the maintenance of 
the V • B = constraint. It should be noted that a solu- 
tion to the latter problem does not necessarily resolve the 
former, since maintaining V ■ B = in a particular numer- 
ical discretisation does not guarantee that it is zero in all 
discretisations. 

W ith regards to the first issue, Bra ckbill fc Barnes! 
(1980) first noted that, when using a conservative for- 
mulation of the magnetic force, a supposed steady state 
could become polluted because of the small but non- 
zero component of magnetic force directed along the field 
lines due to a non-zero V ■ B. This error can have se- 
rious consequences even though the proportional error in 
the magnetic field is small. In SPMHD the force par- 
allel to the field in conservative formulations can have 
catastrophic consequences, le ading to numerical instabi l- 
ity under some circu mstances JPhillips fc MonaghanllT98l . 
iBrackbill fc Barnesl dl98Cfl approached this problem by pre- 
ferring a non-conservative formulation of the momentum 
equation which guarantees that the magnetic force is exactly 
perpendicular to the field. Such an approach has also been 
used successfully in an SPMHD context by several authors 
fe.g.lBenz 1984 ■ Mcdick i et aT]ll995t Efr develd fc Pongrackl 
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merical simulations of shocks seem to require the exact con- 
servation of momentum in order to provide the correct jump 
conditions at shock fronts (which means, at the very least, 
the discrete formulation should be based on continuum equa- 
tions which conserve momentum exactly even with a non- 
zero magnetic divergence). 

This issue of neglect or inclusion of terms proportional 
to V ■ B was discussed at some len gth in paper I , where we 
followed both I Janhunenl ( 2000i) and lDellarl (2001) in formu- 



lating the MHD equations such that they form a consistent 
set in the presence of magnetic monopole terms, retaining 
both the conservation of momentum and energy necessary 
for shocks but using a 'non-conservative' formulation of the 
induction equation. The SPMHD equation set used in pa- 
pers I was shown to form a consistent set with respect to the 
monopole terms in both discrete and continuum forms by 
deriving the SPMHD equations of motion and energy from 
a variational principle which uses the discrete formulations 
of the continuity and induction equations as constraints (see 
paper II). The implications of the formulation of the MHD 
equations in the propagation of divergence errors is discussed 
further in Htj.ll and examined numerically in H7.2I 

Many approaches to the second issue (namely the main- 
tenance of the V ■ B = constraint) are possible. Perhaps 
the simplest in an MHD context is to explicitly evolve a 
vector potential A, from which the magnetic field is derived 
by taking the curl, guaranteeing that the divergence is zero. 
The major disadvantage of this approach is that the com- 
putation of the force terms involves second derivatives of 
the evolved variable (A), which in general can be signifi- 
cantly less accurate. Furthermore evaluation of dissipative 
terms proportional to V 2 B would require computation of 



the third derivatives. Whilst it may be possible to use the 
vector potential in an SPMHD context without degrading 
the accuracy substantially, we do not pursue such an inves- 
tigation in this paper (although it is our intention to do so 
else where) . 

IBrackbill fc Barnesl (Il980l) proposed a simple projection 
scheme to 'clean up' the magnetic field at each timestep, an 
approach which i s now commo nly used in many grid-based 
MHD codes (e.g. Balsara 1998). Similar schemes have been 
implemented in an SPH context for the simu lation of in- 
compressible flows dCummins fc Rudma n 1999). The disad- 
vantage of this approach is that it involves the solution of 
a Poisson equation which is computationally expensive. In 
self-gravitating SPMHD this approach holds some promise 
as the cost may be mitigated by utilising the treecode used 
in the calculation of the gravitational force. In this paper we 
examine various projection methods based on this approach 
in EH 

Another approach used in grid-based MHD codes is 
the so-called 'constrain ed transport' method pioneered by 
lEvans fc Hawlevl 1^)88) in which differences of the mag- 
netic field across the grid cell are constructed in such a way 
as to maintain the divergence free condition exactly. Such 
methods work very well, but is difficult to see how they 
can be made applicable to SPH because of the absence of 
a spatial grid (although perhaps some divergence-free inter- 
polation could be devised). A comparison between several 
constrained -transport type sch emes with the source term 
approach of lPowell et alj <|l999l) and the pr ojection method 
has been recently presented bv lTothl (|2000) for finite differ- 
ence codes. Although not all of the schemes are applicable 
in an SPH context, many of t he nu merical tests presented 
in this chapter are taken fromlTothl's pa per. 

More recently iDedner et alJ (l2002f) have proposed a 
method for cleaning the magnetic field which is significantly 
faster than the projection method. This method involves 
explicitly adding a constraint propagation equation which 
is coupled to the evolution equation for the magnetic field. 
This equation propagates the divergence error in a hyper- 
bolic (ie. wave-like) manner away from its source, allowing 
diffusion of the error to proceed rapidly within the timestep 
condition. This method is easily applied to the SPMHD al- 
gorithm and we provide details of the implementation in 



The paper is organised as follows: In Sj^Jand SjUJwe sum- 
marise the formulation of the continuum MHD equations 
and the corresponding SPMHD form of these equations from 
papers I and II. In the course of the multidimensional test- 
ing, several aspects of the algorithm have been changed or 
refined from that presented in papers I and II. The first 
change is the method for removing the tensile instability, 
which is therefore discussed in Jl| The implementation of 
the dissipative terms formulated in paper I in order to cap- 
ture shocks (paper I) is reviewed and modified accordingly in 
S|K] In Sjjwe investigate several of the approaches discussed 
above to the maintainance of the V ■ B = constraint which 
are applicable in an SPH context, namely the source term 
approach discussed in the p revious chapte r ( H6.1H . projec- 
tion methods f H6.2t and the lDedner et all approach ( H6.3H . 
The algorithm is then benchmarked, as in the one dimen- 
sional gainst a wide range of multidimensional test 
problems used to test recent grid-based MHD codes (aZJl. 
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The tests involve the propagation of an initially non-zero 
magnetic divergence ( H7.2L nonlinear Alfven waves ( H7.3H . 
two dimensional shock tubes ( WJ.b\ . an MHD rotor problem 
(EHJ an d the Orszag-Tang vortex ( H7.7I . The results are 
summarised in 3H1 



The equation set is closed by an appropriate equation 
of state, which for an ideal gas is given by 



(7 - l)pu, 



(9) 



where P is the pressure, u represents the internal energy per 
unit mass and 7 is the ratio of specific heats. 



2 CONTINUUM EQUATIONS 

Our SPMHD formalism is based on the equations of Mag- 
netohydrodynamics in the form 
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s a = _ PS a + _L ( B i B o _ i d lJ B 2 

pa V 2 

The formulation of these equations with respect to terms 
proportional to the divergence of the magnetic field was dis- 
cussed in paper I and derived self-consistently from a varia- 
tional principle in paper II. The implications of this particu- 
lar formulation of the continuum equations in the propaga- 
tion of divergence errors is discussed in HbMl and confirmed 
in the numerical tests presented in 

Note that in place of the specific total energy e, the 
thermal energy can be evolved according to 



du 
It 



P dv i 
p dx* 



(7) 



Similarly, in place of (J1J we could equivalently evolve the 
magnetic flux density B 1 according to 



dB 



= B > dvJ , B i dv \ 

dt dxi dxi 



The difference between evolving the total energy e in 
place of the thermal energy u is found to be very minor. 
One disadvantage of using the total energy is that it does 
not guarantee a positive thermal energy (although this can 
be a useful diagnostic of when a simulation is going wrong) . 
We choose to evolve the magnetic flux per unit mass B l / p 
since it is the natural variable to be carried by particles 
of fixed mass. Again, however, the difference between using 
Q and is found to be minor (although some difference 
might be expected for simulations involving large changes in 
the density) . In general the differences between evolving dif- 
ferent variables in SPH (and SPMHD) is dependent purely 
on the timestepping algorithm used and can be shown to 
decrease as smaller timesteps are used. It should be noted 
that these differences are much smaller than those found in 
grid-based codes due to the exact treatment of the advection 
terms in Lagrangian formulations. 



3 SPMHD EQUATIONS 

The discrete formulation of the SPMHD equations was dis- 
cussed in paper I and derived self-consistently from a vari- 
ational principle in paper II. A self-consistent formulation 
of the SPMHD equations in the case of a variable smooth- 
ing length was also derived in paper II. We summarise the 
equations describing the method below. 

The continuity equation is expressed by the density 
summation 



pa = 2^ m bW(\r a - r b \, h a ), 



(10) 



where W(jr a — ^b\,h) is the interpolation kernel with 
smoothin g length fe ; for w hich we use the usual cubic spline 
(paper I. lMonag han 1992). The time derivative of (IlUt gives 
the SPH expression for the continuity equation (0, in the 
form 



dp a 
dt 



(11) 



where w a b = v„ — vj, and 51 is a normalisation term which 
takes account of the variation of the smoothing length with 
density (paper II), given by 



Ma 



dh a 
dp a 



dWa b (h a ) 



dh a 



(12) 



The smoothing length is assumed to depend on the density 
via the relation 



h a = T) 



Pa 



1/u 



with derivative 
dh a _ h a 

1 ; dp a Vp a ' 



(13) 



(14) 



where v is the number of spatial dimensions. We enforce 
this relation by calculating both h and p self-consistently 
by iteration of the density summation 1101 . The manner by 
which this is done is d escribed in more detail in paper II 
and also in iPricd i2004f) . In brief, since the density of each 
particle is independent of neighbouring particle smoothing 
lengths we are able to iterate only those particles which have 
not converged. Furthermore we predict the new value of the 
smoothing length using the time derivative 



dh a _ h a dp 
dt vp a dt 

As a result the additional cost involved is minimal. 



(15) 



1 Beware that the expression given for C in paper II contains 
some typographical errors. The correct expression is given by 1121 . 



4 Price 



The momentum equation J5J in SPH form is given by 
S lJ \ dW ab (h a ) 

np 2 ) a dxi 



if = 2>» 



S lJ \ dW ab (h b ) 



Vp 2 J b dxi 



(16) 



whilst the energy equation in discrete form is given by 



dCa 

dt 



= E 



nn, 



tip 2 



v b V 3 a W ab {h a ) 
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n P 2 



(17) 



The internal energy equation @ in SPH form is given by 

(18) 



Finally, the induction equation is given by 

d_( B i 
dt 



B j i dW ab (h a ) 

) _ m b v ab - 



dxi 



or alternatively 
dt Q, a p, 



— V m b \v l ab B 3 a - B l a v 3 al 

a Pa , L 



dW ab (h a ) 
dxi 



(19) 



(20) 



4 INSTABILITY CORRECTION 

In paper I an artific i al str ess or 'antidumping' term de- 
scribed bv lMonaghanl jl997l) was used to eliminate the ten- 
sile instability associated with a conservative formulation of 
the momentum equation in SPMHD. The basis of this ap- 
proach is that the instability manifests as particles clump- 
ing together in the presence of a negative stress and at 
short wavelength s (see below). The solution proposed by 
iMonagharJ il997f) and described in paper I was to introduce 
a repulsive term proportional to the anisotropic magnetic 
force which acts to remove the instability at short wave- 
lengths by preventing particles from clumping together. This 
term has been used very e ffectively in elastic dynamics sim- 
ulations llGrav et alJl20o3) and was found to remove the in- 
stability very effectively in the one dimensional simulations 
considered in papers I and II. A more detailed investi gation 
of the antidumping term has been given recently in iPricd 
( 2004), interpreting the antidumping term as a modification 
of the kernel gradient used in the anisotropic force term. In 
this investigation several disadvantages to the antidumping 
approach were highlighted. The first is that at large negative 
stresses (e.g. at low magnetic /3) the antidumping term can 
cause the numerical estimate of the sound speed to be sig- 
nificantly in error. The second, somewhat fatal disadvantage 
is that the antidumping term does not appear to guarantee 
stability in the case of a variab le smoothing length. For more 
details we refer the reader to iPricel <l2004h . however it suf- 
fices to say that the antidumping approach does not appear 
to be uniformly satisfactory for dynamical MHD problems, 
particularly in more than one dimension. We therefore con- 
sider two alternative approaches in this paper, which are 
outlined below ( H4. 114.21 . 

It is worth recalling that the physical source of the in- 
stability is the additional small but non-zero force directed 



parallel to the magnetic field in the conservative formulation 
of the MHD equations (see introduction) . For this reason it 
might be expected that enforcing the V ■ B = condition 
might also eliminate the tensile instability. In fact this is 
not the case, since the instability manifests even in one di- 
mension (where the divergence is zero exactly). The reason 
for this is that although the divergence is zero by virtue of 
B x — const, the gradient of this constant (as evaluated in 
the force term) is not necessarily zero numerically. In par- 
ticular this is the case in the conservative formulation of 
the SPMHD equations, since the symmetric SPH gradient 
evaluation in the form 



. V Pa P b 



vw ab 



(21) 



is non-zero in the case of A = const. To counter this problem 
two approaches may be taken. The first is to add or subtract 
an arbitrary constant in order to keep the total stress pos- 
itive and thus preventing negative stresses (and instability) 
from occurring. The second approach is to use an SPH gra- 
dient operator which vanishes for constant functions. These 
approaches are described below. 



4.1 Removing the constant component of 
magnetic field 

A simple method for removing the tensile instability is to 
subtract an arbitrary constant from the stress in order to 
make the total stress positive. For simulations where the 
magnetic field is strong due to an initial net flux through 
the system, a natural choice for this constant is to subtract 
the external (ie. produced by currents outside the simulation 
domain) component of the magnetic field. In this case the 
stress tensor @ for particle a is modified according to 



Pa + ^-Bl) 5 ij + — ( BlBl - BqBq 
2^o / 



(22) 



where Bo is the magnetic field component which does not 
change throughout the simulation (for example in one di- 
mensional simulations we would use Bo = [B x ,0, 0]). In 
general the constant field could also have a spatial profile 
(for example in a fixed dipole field from the central star in 
an accretion disc) in which case the analytic gradient could 
be used. In all of the cases we consider the external mag- 
netic field is always the same independent of the particle 
position, such that calculating l|22|l involves storing only a 
single vector. It is worth noting that the formalism given 
above (where the constant field is subtracted from the total 
field) is more efficient than explicitly adding the contribu- 
tions from separate constant and variable field components. 

This simple solution completely cures the one dimen- 
sional instability because the B x component of the field is 
explicitly removed from the anisotropic gradient term. Neg- 
ative stresses can only arise in this formulation when the 
anisotropic terms in the fluctuating component dominate 
the isotropic pressure term (from which the constant field 
has not been subtracted). 

A more general formulation which guarantees stability 
at all times must ensure that the total stress is positive. This 
can be achieved by using 



5,'; 



Scons 



(23) 
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where 

Sconst — max 



2~p7^ 



.0 



(24) 



where the maximum is taken over all the particles. In 
many ways this is s imilar to the original proposal of 
iPhillips fc Monaehanl il985T> in which the maximum value 
of the stress tensor over all the particles was determined and 
then subtracted from the stress for each particle. One disad- 
vantage to this approach is that total energy is not conserved 
exactly since the contribution to the total energy evolution 
from the induction equation (which uses the total magnetic 
field) does not exactly balance the contribution from the mo- 
mentum equation. An alternative is the approach of Morris 
(described below) in which the anisotropic term is modi- 
fied slightly. In this paper we subtract the external field in 
simulations where a dominant external field is present, as in 
many of the two dimensional problems considered in this pa- 
per, reverting to the Morris approach otherwise. In practice 
we find little to differentiate the two approaches. The results 
in all cases are much better than those obtained using the 
antidumping term. 



4.2 Morris approach 

An approach suggested by iMorrisI il996T) is to retain the 
conservation of momentum on the isotropic terms in 11161 but 
to treat the anisotropic terms using a differencing formalism 
which is exact in the case of a constant function (see above). 
The force term is then given by 



E 



nib 



Pa + \Bll^ P b + \BH^ \ OWa. 



pi 



+ 



-E 



m b - 



(BjBj)b — (BiBj) a dWgb 
paPb dXj 



(25) 



(26) 



This formalism does not therefore guarantee exact momen- 
tum conservation (since the anisotropic term does not give 
equal and opposite forces between particle pairs) but can 
be expected to give good results on shocks for which the 
anisotropic term is less important. It is also a better ap- 
proach than using formalisms based on a pure J x B force 
since (1261 is still a discretisation of a tensor force and there- 
fore conserves momentum in the continuum limit for non- 
zero V ■ B. This also means that 1261 retains the consis- 
tent formulation of the MHD equations in the presence of 
monopoles, although the discrete equations are no longer 
self-consistent with each other (where self-consistent means 
that the equations can be derived from a variational princi- 
ple and will thus conserve momentum, energy and entropy) . 
Note that when using the variable smoothing length terms, 
we use the average of the normalised kernel gradient in 
12611 . as in the dissipative terms. The small amount of non- 
conservation introduced by the Morris formulation is not 
found to significantly affect the shock capturing ability of 
the scheme. 



5 DISSIPATIVE TERMS 

Artificial dissipation terms which are required in order to 
simulate shocks were formulated in paper I. These terms 



are given by 



dv\ 
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dB a 

dt 

de a 
dt , , 
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A _ (r n - r b ) 
|r a - r b \ 



2_^m b VaWab, (27) 



Pab 



= P 
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a ^2 m b aB _2 Sl3 (r x (B ab x r)r ab F ab {28) 

b 

-E m6 " 



Pab 



r • VaWab, (29) 



(30) 



and 



|a(v a • f) 2 + . 



e *a = { +\ a B[Bl - (B a • i) 2 ]/nop ab 



Vab ■ Tab < 0; 



a u u a + I ct B [Bl - (B a ■ rf]/p Pab, v ab -r ao >0; 

with a similar expression for el . Note that the notation used 
in this paper differs slightly from that used in paper I. In 
particular we use separate parameters a, ct u and as to con- 
trol the artificial viscosity, thermal conductivity and resis- 
tivity respectively rather than a single parameter K. Note 
also that these parameters are expected to be of order unity 
rathe r than K ~ 0.5 (su ch that a corresponds to the a used 
in the lMonaghanl il992Tl artificial viscosity formulation used 
widely in SPH). 

The signal velocity v 3 i g represents the fastest speed of 
signal propagation between the two particles. In MHD we 



[Va +Vb~ /3v ab ■ f ] 



(32) 



where 



Va = ~ 



cl + 



B 2 

^a 
PapO 



+ 



2B Q 



VPaPO 



+ 



2B a ■ rc 



B 2 



(33) 



with a similar equation for Vb, where c is the sound speed. 
Again our notation differs slightly from that used in pa- 
per I as we use a dissipation parameter a of order unity and 
Vsig ~ c s (as opposed to K ~ 0.5 and v S i g ~ 2c s in paper I). 
The ft term in the signal velocity in this formalism naturally 
provides the non-linear (Von Neumann-Richtmyer) compo- 
nent of the artificial viscosity. 

The dissipative terms 1271 and 1281 provide an artificial 
viscosity and resistivity. The term involving (u a — U b ) in the 
energy equation provides an artificial thermal conductivity. 
These terms are derived so as to guarantee a positive definite 
contribution to the thermal energy and entropy (paper I). 
For reference the dissipative term added to the internal en- 
ergy equation is given by 



E 



m b - 



pab 



-a [(v a • f) - (vf, • r)] 



+a u (u a — u b ) 
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' Bib- (B ab -v) 2 } 



2 pa Pa 
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In the one- dimensional tests described in paper I, it was 
found that the dissipative terms as described above were in- 
sufficient in shock problems involving jumps in the trans- 
verse velocity component, resulting in a modification of the 
viscosity term in order to apply the viscosity to all velocity 
components. This was attributed to the fact that the prob- 
lems considered involved two and three dimensional velocity 
components whilst restricting the particles to move in only 
one spatial dimension. In the course of the multidimensional 
tests, it became apparent that this conclusion was incorrect. 
In paper I, following the usual procedure for hydrodynamical 
SPH, the dissipative terms were applied only for approach- 
ing particles (v a (, ■ f < 0). Whilst this is appropriate for 
the artificial viscosity term, discontinuities in the magnetic 
field (requiring artificial resistivity) can occur for particles 
in both compression and rarefaction. Applying the artificial 
resistivity 12811 uniformly across the simulation was found to 
correct the oscillations observed in the magnetic field, and 
hence also in the transverse velocity components, removing 
the need for any modification of the viscosity term. In fact 
the results with the artificial resistivity term applied sepa- 
rately to the viscosity are an improvemen t on those given in 
paper I and are presented in lPricei i2004h . 



5.1 Dissipation terms using total energy 

A further issue in a multidimensional context is that in the 
derivation of the above dissipative terms (paper I) it was 
assumed that only components of the magnetic field per- 
pendicular to the line joining the particles would change 
at a shock front. However, in a multidimensional simula- 
tion the assumption of non-zero magnetic divergence may 
not hold exactly, as has already been discussed. In partic- 
ular divergence errors are often created at flow discontinu- 
ities where fluid quantities are changing rapidly. It therefore 
makes good sense to drop the assumption of non-zero mag- 
netic divergence in the derivation of the dissipative terms. 
The assumption that only the velocity components paral- 
lel to the line joining the particles will change is also not 
strictly true in MHD since velocity components transverse 
to this line will change with a jump in the transverse mag- 
netic field. For this reason we re-derive the dissipative terms 
with an energy term of the form 



1 2 B z a 
-orv a + a u u a + «b- — — 

2 ZfloPab 



(35) 



which involves both the total kinetic and magnetic ener- 
gies. The implication is that smoothing is then also applied 
to jumps in B parallel to the shock (ie. V ■ B jumps) and 
via the kinetic term to transverse velocity jumps (ie. shear 
discontinuities). For the contribution to the entropy to be 
positive definite, the terms in the thermal energy equation 
must take the form 



OB 



pab 



■(B a ~B b ) 2 



2floPab 

+ a u (u a - U b )} r ■ VaWab, 



which correspondingly requires dissipation terms in the mo- 
mentum and induction equations of the form 



dv a 
dt 

dB 

dt 



t-^ av s ig(v a - Vj,)„ 

2_^m b r-VaWab, (37) 

. Pab 



Pa m b ^^ (B„ - B 6 ) f ■ V a H/ ai (38) 

h Pab 



In the multidimensional case we find that use of 1381 has 
distinct advantages over 12811 since in more than one dimen- 
sion divergence errors can cause the extra component of the 
magnetic field to jump slightly. Whether or not to use 137> 
in place of 1271 is slightly less clear. The application of dis- 
sipative terms to spec ific discontin uities in a hydrodynamic 
context is discussed in lPricd (|2004) with regards to artificial 
thermal conductivity, where it was found that smoothing 
of discontinuities in the thermal energy was necessary only 
where the discontinuity is not already smoothed by the ap- 
plication of artificial viscosity (which could occur, for exam- 
ple at a contact discontinuity). In the present case, since a 
jump in transverse velocity can only occur at a correspond- 
ing jump in the transverse magnetic field, these discontinu- 
ities will already be smoothed by the application of artificial 
resistivity there and so the use of 13711 may simply result in 
excessive dissipation (since it must also be applied to parti- 
cles in both compression and rarefaction, whereas the usual 
viscosity term is applied only to particles in compression). 
Furthermore 1371 no longer conserves angular momentum 
(since the viscosity is not directed along the line joining the 
particles) and also no longer vanishes for rigid body rotation 
(since in effect rotational energy is converted into thermal 
energy). Thus for simulations involving significant amounts 
of shear (for example in accretion discs) the effects of using 
137H would need to be studied quite car efully. It is w orth 
noting that a similar term was used by iMorrid dl996Tl for 
SPMHD shocks in place of an artificial resistivity. 



5.2 Dissipation switches 

The artificial viscosity para meter a is controlled usin g the 
switch described in paper I IIMorris fc MonagharJll997 

da a — c 
~dt ~ 7 



+ 5 



(39) 



In paper I we effectively used the source term 
S = max(-2V • v,0.0) 



(40) 



wh ich (as noted in pape r I), is double the source term used 
by|M 

orris fc MonagharJ l|1997l) (note that in paper I a dis- 
sipation parameter K is used which has a value of a/2). In 
paper I it was found that the stronger source term was nec- 
essary in order to effectively damp post-shock oscillations. 
However for certain problems this source term could cause 
the artificial viscosity to become overly strong, resulting in 
excess smoothing of shock fronts. For this reason, in this 
paper we adop t a m odification of the switch proposed by 
iRosswog et alJ (1200(1) , where the source term is given by 



max(-V ■ v,0.0)(2.0 - a) 



(41) 



(36) 



resulting in an initially stronger source term which tails off 
as a reaches its desired value of unity at the shock. 
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Since artificial resistivity is required at discontinuities in 
the magnetic field, which may occur where particles are not 
necessarily approaching each other, artificial viscosity and 
resistivity should not be controlled using the same switch, as 
was the case in paper I, leading to unnecessary modifications 
of the artificial viscosity term. A similar switch appropriate 
to the artificial resistivity term can be devised similar to 
that used in the viscosity. We evolve the resistive dissipation 
parameter ob according to 

dots 
dt 



CtB 



+ s 



(42) 



where in this case the source term is given by 



5 = max 



(43) 



which has dimensions of inverse time, as required by 1421 1. 

A similar switch may also be derived for the artificial 
thermal conduct ivity. A swit ch based on the first derivative 
of u was used in IPricd <l2004f) . In this paper we use a switch 
based on the second derivative of it, where the source term 
is given by 



5 = 0.1/ilV u\ 



where h is the smoothing length and we multiply the source 
term by a small number in order to apply only the very min- 
imum amount of dissipation needed to eliminate the wall 
heating effect. The second deriva tive term is computed ac- 
cording to (e.g. lBrookshawlll985l) 



- U b ) r ab ■ VgWa 

Pb r 2 ab 



(45) 



The second derivative switch is preferable since it responds 
only to sharp discontinuities in u, ensuring that a minimal 
amount of artificial thermal conductivity is applied. Ad- 
ditionally it requires storage of fewer quantities than the 
switch involving the first derivative. We have not investi- 
gated the use of switches for artificial viscosity or resistiv- 
ity based on second (or higher) derivatives in this paper 
although it deserves further study. In particular a switch 
based on V(V • v) would be very useful in self-gravitating 
simulations where —V • v can have large constant values in 
the absence of shocks due to the gravitational collapse. 



6 DIVERGENCE CORRECTION 
TECHNIQUES 

6.1 Source term approach 

The induction equation can be written in the 'conservative' 
form 



f = -Vx(vxB), 



(46) 

= V-(vB-Bv). (47) 
which explicitly conserves the volume integral of the flux 

BdV (48) 



In Lagrangian form 1461 can be written as 
= _ B (V ■ v) + (B ■ V)v + v(V ■ B). 



(49) 



Taking the divergence of this equation, we have 
B) =0, 



(50) 



showing that the constraint V ■ B = enters the MHD 
equations as an initial condition. However allowing mag- 
netic monopoles resulting from V • B / to evolve ap- 
propriately within the flow can prevent the build up of un- 
physical numerical effects associated with their presence and 
can therefore reduce the need for comp utationally ex pen- 
sive divergence cle aning procedures. Thus lPowelll (1 19941) (see 
iPowell et allll999l) suggested that the conservative forms of 
the MHD equations should contain source terms to ensure 
that the se errors are p ropagated out by the flow. With this 
in mind, IPowelll Jl994T) added source terms to the momen- 
tum, energy and induction equations, which take the (La- 
grangian) form 



1 dS l - 



B i dB J 



<M_ 

dt p dxi p dx 3 ' 

de 1 d(viS i j ) v t B l dB 3 

dt p dx 3 



dx 3 ' 



dv 3 



(44) -7T - BJ »ZJ- B 7^- 



(51) 
(52) 
(53) 



dB_ 

dt " dxi " dxi 

Taking the divergence of (1531 shows that the divergence er- 
rors in this formalism evolve according to 



^(V-B)-rV-(vV-B) 



0, 



(54) 



which has the same form as the continuity equation for the 
density (where in this case we have a density of magnetic 
monopoles, V ■ B). This therefore implies that the total vol- 
ume integral of V • B across the simulation is conserved and 
hence that the surface integral of the flux 



B ■ dS = / (V ■ B)dV, 



(55) 



is conserved. The conservation of this quantity is far more 
important physically than the conservation of the volume 
integral (Blfl l. 

The disadvantage of using I5H - 153| I is that exact conser- 
vation of momentum and energy is sacrificed, which proves 
to be important for shock-type problems. Correspondi ngly it 
can l ead to incorrect jump conditions at s hock fronts jT^thJ 
l2000jl. More recen tly it has been shown bv ljanhunenl (2000) 
and [Dollar (2001) that the correct formulation of the MHD 
equations in the presence of monopoles should not violate 
the conservation of momentum and energy. 

The 'mo nopole formulation' of Ijanhunenl ll2000f) and 
Dollar i200ll) is identical to the self-consistent formulation 
of the SPMHD equations derived in paper II and given by 
Q-(0J|. Note that the induction equation JSJ (equivalently 
using Q) is the same as in Powell's method and therefore 
the same conclusions can be drawn regarding the manner in 
which the divergence errors evolve (1541 . We investigate the 
implications of the 'source terms' in the induction equation 
in 

6.2 Projection methods 

A common approach to the divergence problem is to clean 
up the magnetic field at regular intervals via the projec- 
tion method fe.g. iBrackbill fc Barneslll980l) . The basic idea 
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is to decompose the magnetic field into a curl and a gradi- 
ent (which can be done unambiguously for any vector field) 
according to 



B* = V x A + V0. 



(56) 



From this decomposition there are two ways of obtaining a 
divergence free field, both of which we discuss below. 



6.2.1 Scalar projection 

Taking the divergence of this expression results in the Pois- 
son equation 



V 2 </> = V-B*, 



(57) 



which can then be solved for the scalar quantity <j>. The 
magnetic field is then corrected according to 



B = B* - V<j>. 



(58) 



The major disadvantage with this approach is that the solu- 
tion of the Poisson equation 157H is computationally expen- 
sive, scaling as 0{N 2 ). In an astrophysical SPH context this 
may be offset somewhat by the fact that the Poisson equa- 
tion for t he gravitational field is u s ually solved using a tree 
code (e.g. iHernauist fc Katdllflsll iBenz et al.lll99(t l which 
scales as O(iVlogiV). There are, however, some subtleties to 
this approach, which we outline below. 

Projection schem es for incompressible flow in S PH have 
been implemented bv lCummins fc RudmarJ il999Tl . the re- 
sults of which are applicable to the present case. The im- 
portant point, also discussed bv lTothl (2000) is that for the 
projection step to reduce the divergence to zero (ie. to pro- 
vide an exact projection) requires that the discrete version 
of 15711 is satisfied exactly. This means that the operator 
used to evaluate the divergence term on the right hand side 
of 15711 should be the same as the divergence operator used 
in the evaluation of the V 2 on the left hand side and that 
the gradient operator used in t he evaluation of V 2 should be 
the same as that used in 158t . ICummins fc RudmarJ (ll99St) 
approach this problem by calculating the V 2 using SPH op- 
erators, solving the Poisson equation by matrix inversion. 
Good results were also obtained using an approximate pro- 
jection (ie. where the divergence ope rators on the left and 
right h and side differ). In this scheme ICummins fc RudmarJ 
( 1999) used the SPH evaluation of the Laplacian similar to 
that which is commonly used for thermal cond uction in SPH 
jBrookshawlll985l:IClearv fc Monaghanlll999l) (and which is 
similar to the artificial thermal conductivity term used in 
this paper). The Poisson equation is then solved by invert- 
ing the resulting matrix equation. 

The solution of 15711 by direct summation (of which the 
tree code is an approximation), uses the exact solution to 
the Poisson equation 1571 given by 

4>{t) = j G(|r-r'|)V • B(r')dV(r'), (59) 
where G(|r — r'|) is the Green's function, given by 
G(r) — — In r + const, 

°W = -4^' (60) 



in two and three dimensions respectively. The gradient 
needed in the correction step can be calculated directly, giv- 
ing (in three dimensions) 



V*(r) = -£ 



V • B(r') , , , 
— — — (r-r )dV(r ). 



(61) 



In SPH we replace the volume element pdV with the mass 
per SPH particle and write the integral as a summation 
according to 



V(^ a = — 22 1Tlb 

b 



(V ■ B) b (r a - r b ) 

4:Trp b \r a — r b | s 



(62) 



Since we still retain the freedom to choose the discrete oper- 
ator used to evaluate V • B at each particle, it becomes clear 
that the solution by direct summation will only provide an 
approximate projection, since l)57|l is not discretely satisfied. 
This approximate solution will be degraded further when im- 
plemented using a tree code. A further disadvantage of the 
projection method for many of the problems considered in 
this paper is that it is somewhat complicated to implement 
in the case of periodic boundary conditions. The implication 
of these subtleties in the practical application of the projec- 
tion method based on the Green's function solution (using 
a direct summation over the particles) are discussed in H7.2I 
Essentially we find that this projection method is reason- 
ably effective at removing divergence errors at wavelengths 
larger than the smoothing length, but is less effective at re- 
moving short wavelength (~ h) noise due to the smoothing 
of this noise inherent in the SPH operator used to calculate 
V-B. Preliminary calculations using this projection method 
in conjunction with a tree code in three dimensions indicate 
similar results. In this paper we compute the divergence of 
the magnetic field using the SPH operator 



(V-B) 



1 



J2 m b (B a - Bfc) ■ V a Wab(h a 



(63) 



6.2.2 Vector projection 

An alternative projection scheme can be implemented by 
solving for the vector potential A. That is, we take the curl 
of 156H to obtain 



V x B* = V(V ■ A) — V 2 A. 



(64) 



Choosing the Gauge condition V- A = 0, we obtain a Poisson 
equation for the vector potential in terms of the current 
density J = V x B*//xo 

V 2 A = -At J (65) 
with solution 



A(r) 



G{\r- r'|)J(r')dV(r'). 



(66) 



Taking the curl, we obtain an equation for the corrected 
magnetic field in terms of the current density, which in three 
dimensions is given by 



47T 



J(r') x (r-r') 



dV(r') 



(67) 



which is simply Biot-Savart's Law. In SPH form this is given 
by 
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B a = -£ 



(V x B*) h x (r a - r b ) 
4"Kp b \r a -n\ 3 



(68) 



This method could also be useful in an SPH context in sit- 
uations where several disconnected regions exist contain- 
ing strong magnetic currents. By solving 16 7H . the cor- 
rected magnetic field is determined from the current den- 
sity, resulting in a knowledge of the magnetic field at any 
point in space. This approach was in fact used as the ba- 
sis for the very first SPMH D algorithm implemented by 
iGingold fc Monaghanl jl977l) . In this paper we will only con- 
sider the use of 1681 as a divergence cleaning method. In this 
respect solving I65H has a slightly higher computational ex- 
pense than 1621 since the Poisson equation is solved for a 
vector quantity rather than a scalar, giving (up to) three 
summations of 0(N 2 ) as opposed to just one. Nevertheless, 
there is a significant difference between the two methods. 
The difference is that whereas the approximate nature of 
the projection in 16211 means that the divergence is not guar- 
anteed to be reduced to zero, in 16811 the divergence of the 
expression for B is zero exactly by virtue of the curl in the 
summation. 

The approximate nature of the projection in this case 
means that the current is not guaranteed to remain un- 
changed during the projection step. However (as noted by 
lMonagharJl992h the current is usually well estimated by the 
SPH particles since the current is in general where the mat- 
ter is. In this paper we compute the current density using 
the SPH operator 



(V x B)„ 



1 



J2m b (B a -B b ) xV a W ab (h a ). (69) 



In practise we find that this projection method is far 
more effective than the scalar projection and this is demon- 
strated in H7.2I fsee Figure |5J. Again preliminary three di- 
mensional calculations indicate similar results, although an 
implementation using the tree code is more difficult in this 
case since all three components of the vector potential must 
be stored and summed over the tree as opposed to just one 
in the scalar projection (in which case the standard grav- 
ity tree can simply be called with the source term replacing 
the particle mass). However the degree to which the physi- 
cal current is affected by this projection in three dimensions 
remains to be investigated. 



6.3 Hyperbolic divergence cleaning 

iDedner et alJ (|2Q02) examine alternativ e divergence clean - 
ing procedures. In their paper (see also lMunz et, al.l feoOO'l. 
they derive a general constrained formulation of the MHD 
equations, from which formalisms can be derived to give di- 
vergence cleaning which is elliptic (involving the solution of 
a Poisson equation), parabolic (in which the divergence er- 
rors are diffused away) and hyperbolic (where the divergence 
errors are propagated away from their source at a character- 
istic speed). The projection method described above is an 
elliptic approach, the main disadvantage to which is the sub- 
stantial computational cost involved in the solution of the 
Poisson equation. The parabolic approach was found to be 
severely limited in scope due to the timestep restrictions im- 



posed by the Courant condition 2 . The hyperbolic approach 
was found to be particularly effective, especially when com- 
bined with a parabolic term such that divergence errors are 
both transported and diffused. It is this approach that we 
outline below in an SPH context. 

The basic idea is to introduce an additional scalar field 
ip, which is coupled to the magnetic field by a gradient term 
in the induction equation, 

^ = -B(V ■ v) + (B ■ V)v - Vnji. (70) 
at 

Note that our induction equation maintains the consistent 
treatment of divergence terms discussed above. The variable 
ip is then calculated by adding an additional constraint equa- 
tion, which for the combined hyperbolic/parabolic approach 
is given by 

*_-4<v.b,-4. m) 

Neglecting the second term on the right hand side of (1711 
gives an equation for ip which is purely hyperbolic. This im- 
plies that divergence errors are propagated in a wave-like 
manner away from their source with char acteristic speed Ch 
(for more details we refer the reader to the lDedner et all pa- 
per). The second term on the right hand side is a parabolic 
term which causes ip to decay exponentially to zero with 
c-folding time r (this is easily seen by neglecting the hy- 
perbolic term and solving the resulting ordinary differential 
equation for ip(t)). Since it is desirable for the divergence er- 
rors to be propagated at the maximum possible rate (within 
the timestep constraint imposed by the Courant condition), 
Ch should be set equal to the maximum signal propagation 
speed. For simplicity we calculate this as 



Ch = 



7 P 1 B 2 
p 2 pop 



(72) 



where the maximum value over all of the particles is used. 
The gradient term in the induction equation is calculated 
using a simple SPH estimate 



Qa.Pt 



22 m htyb ~ V'a)Va W ab (h a ) 



(73) 



Similarly the divergence of the magnetic field is calculated 
using I63H . 

The choice of d ecay timescale r is more complicated. In 
IDedner et alJ i2002f) the decay timescale used is given by 



Ch 

c r 



(74) 



where they find that an optimal cleaning on their chosen test 
problem is given by c r — 0.1. The problem with this is that 
c r is not a dimensionless parameter, but rather has units of 
length. Thus the optimal choice for any given problem will 
depend on the length scales in that particular problem. We 
explicitly write the timescale as 



1 

Ta 



(JCh 
Xa ' 



(75) 



2 an equivalent approach in SPMHD is to use an artificial resis- 
tivity in order to diffuse away divergence errors. This has been 



used, for example, bv lMorrisI il996l) and lHoskingl i200 
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where A is a length scale and a is a dimensionless param- 
eter which determines the decay timescale. Setting a = 
therefore gives a purely hyperbolic correction. The physi- 
cal interpretation of the length scale in the problem can be 
determined by solving the following reduced system of equa- 
tions (ie. neglecting the usual MHD evolution terms) 



dB 

dt 



— = -v</> 



— - -4(VB)--. 



(76) 
(77) 



Taking the divergence of the (1761 and substituting into 177H 
we obtain the following equation for ip 



1 d 2 i> 2 aX dip 



(78) 



where an identical equation may be obtained for the evo- 
lution of V • B. This equation is simply the wave equation 
with a damping term, the solution to which is easily obtained 
by a separation of variables and is given in many standard 
textbooks. The length scale enters the solution as the wave- 
length for critical damping, that is the wavelength at which 
solutions change from being wave-like to being damped. 

In practical calculations we expect divergence errors to 
be generated at wavelengths close to the smoothing length. 
We therefore set X — h and determine the value of the di- 
mensionless parameter a by experiment. A value of a — 0.2 
would imply that tjj (and thus V ■ B) will have decayed sig- 
nificantly after the divergence errors have propagated ap- 
proximately 5 smoothing lengths. In H7.2l we examine in de- 
tail the effects the hyperbolic cleaning on a problem involv- 
ing a fixed wavelength of error (ie. independent of h) which 
graphically illustrates the divergence cleaning method (see 
Figure [3J. The effect of this type of cleaning on errors gen- 
erated by the flow are examined in H7.7I We find that values 
of a ~ 0.4 — 0.8 generally give the best results, giving a good 
balance between the hyperbolic (fast but non-diffusive) and 
parabolic (diffusive but slow-acting) effects. In practise some 
diffusion is also added by the artificial resistivity terms (©. 
In general, however, the divergence correction provided by 
the hyperbolic/parabolic scheme is found to be quite small 
(around a factor of ~ 2 reduction). Thus, whilst this type 
of divergence cleaning essentially comes free-of-charge com- 
putationally, under some circumstances it may be necessary 
to supplement it with a stronger form of cleaning, such as 
use of a projection method or some other kind of elliptic or 
parabolic cleaning which is not limited to the explicit time 
step condition. 



7 NUMERICAL TESTS 

The main issue to be addressed in 2D and 3D problems is 
the non-zero divergence of the magnetic field. In the SPH 
context it also allows us to estimate the extent to which the 
artificial dissipation spuriously affects the numerical results. 
Again there is a substantial literature of multi-dimensional 
MHD problems which have been used to test grid-based 
MHD codes (e.g. iDai fc Woodward! 11994 lRyu et alJll995t 
lBalsaralll99l IDai fc Woodward! Il998l : iTotbl l200(tl and we 



7.1 Implementation 

The implementation of the SPMHD equations used for the 
multidimensional tests is almost identical to that used in 
the one dimensional case (paper I). The density is calcu- 
lated by summation, the total energy equation is used (al- 
though results are indistinguishable using the thermal en- 
ergy equation in nearly all cases) and the magnetic field is 
evolved using I19H (or using I70H when using the hyperbolic 
cleaning). In the shock tube tests we use unsmoothed ini- 
tial conditions. The artificial dissipative terms, except where 
otherwise indicated are implemented using the jump in total 
magnetic energy f H5.ll but the viscosity term uses only the 
velocity component along the line joining the particles 12711 . 
Unless otherwise indicated, artificial viscosity and thermal 
conductivity are applied using the switches discussed in H5.2I 
whilst the artificial resistivity term is applied uniformly us- 
ing ob — 1. A major difference between the simulations 
presented here and those in the paper I is that the anti- 
dumping approach was not found to be uniformly successful 
in eliminating the tensile instability for all of the problems 
considered (in particular for the Alfven wave test only a 
narrow range of parameters would produce stable results). 
Furthermore this term was found to result in spurious extra 
numerical noise, particularly in the shock tube tests. For this 
reason we have eliminated the tensile instability by simply 
subtracting the constant component of the magnetic field 
from the gradient term f H4.lt in the shock tube problems 
and using the stable Morris formulation of the magnetic 
force f H4.2H elsewhere. Note that even on the shock prob- 
lems the differences in results between these two methods is 
almost negligible. 



7.1.1 Error estimates 

Various estimates can be made of the error produced in the 
simulation by any non-zero magnetic divergence. Monitoring 
these quantities over the course of a simulation thereby gives 
some measure of the magnitude of the error produced by 
V • B. The most common approach in SPH implementations 
to date has been to monitor the dimensionless quantity 



IB 



(79) 



and ensure that it remains small (typically < 0.01) over most 
of the simulation, where h is the SPH smoothing length and 
the divergence is calculated using 1631 . This provides some 
measure of the relative error in the magnetic field but no in- 
dication of how much influence this error has in the dynam- 
ics. For this reason it is also useful to measure the relative 
error in the total force caused by a non-zero divergence, 



Ef 



force 



IfllBI 



(80) 



consider several of these problems here. 



where f mag is the magnetic component of the SPH force 116H . 
whilst f is the total force on the particle. It is also useful to 
simply monitor the evolution in the maximum, minimum 
and average of | V • B| with time as well as various conserved 
quantities. 
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7.1.2 Conserved quantities 

Aside from the usual conserved quantities of mass, momen- 
tum, angular momentum, energy and centre of mass, several 
additional quantities can be measured in MHD which can be 
useful diagnostics in a numerical simulation. A list of such 
quantities can b e derived using Hamiltonian t echniques and 
is given by fe.g.'l lMorrison fe Hazeltind lll984t) . The helicity, 

y"(A-B)dV, (81) 

where B = V x A, is a measure of the linkage of magnetic 
field lines (expressing the fact that magnetic field lines which 
are initially linked cannot become unlinked in the absence of 
dissipative terms). This quantity can only be usefully mea- 
sured in simulations which explicitly use the vector potential 
A. A similar invariant is the cross helicity 

/(B-vjdV^Vm^-v, (82) 

J V pb 

which measures the mutual linkage of magnetic field and 
vortex lines. The conservation of the cross helicity is a re- 
sult of the magnetic field lines being frozen into the fluid. 
Measurement of the conservation of this quantity in a nu- 
merical simulation therefore provides an estimate of the de- 
gree of slippage of the magnetic field lines through the fluid. 
The volume integral of the magnetic flux 1481 is also con- 
served across the simulation volume, provided that the flux 
is normal to (or zero at) the boundary of the integration 
volume. However the conservation of flux in a volume sens e 
is not particularly important physically dJanhunenl [2000 ) . 
More important is that the surface integral of the flux 1551 
should be conserved. The conservation of these quantities 
with respect to formulations of the MHD equations in the 
presence of magnetic monopoles was discussed in > I6.1I 

There is also a conserved quantity which is the 
MHD analogue of the circu lation (Bck enstein fc Oror]l200ot 
iKuznetsov fc RubanlEoOol) . although the physical interpre- 
tation is somewhat obscure. It has been shown that SPH 
conserves an appro ximate version of the cir culation in the 
hydrodynamic case dMonaghan fc Pricel200lh . related to the 
invariance of the equations to the relabelling of particles 
around a closed loop due to the frozen-in vorticity field. A 
similar, though more restricted relabelling symmetry holds 
in the MHD case (in that the particles around the loop must 
also be on the same field line) and it may therefore be ex- 
pected that SPMHD also maintains this invariance. 

7.1.3 Visualisation 

In order to make a direct comparison of our results with 
those of grid-based MHD codes, we interpolate the results 
from the particles to an array of pixels using the SPH kernel. 
That is, for a contour or rendered plot of a scalar quantity 
<f> we interpolate to the pixels using 

4>{x,y) = V* mj— W(x — Xb,y — Vb, hb) (83) 

where W is the cubic spli ne kernel used in the calculations 
(paper I; lMonaghanlfl992l) and the summation is over con- 
tributing particles. Note that in practise this is quite simple 
to implement, as it involves only one loop over the particles, 



during which the contributions from the current particle to 
all pixels within a smoothing radius (2h) are calculated. For 
a vector quantity a similar interpolation can be performed 
for each component. An interactive plotting program in- 
corporating these interpolation schemes for visualisation of 
SPH (and SPMHD) data in 1, 2 and 3 dimensions has been 
written by the author and is available upon request. 

7.2 V B advection 

The firs t problem we examin e is a simple test similar to that 
used bv lDedner et alJ (120021) in which a non-zero magnetic 
divergence is introduced into the simulation as an initial con- 
dition. This is a particularly good test for comparing var- 
ious divergence cleaning procedures. The initial conditions 
are a uniform density distribution (p = 1) in the domain 
—0.5 < x < 1.5, —0.5 < y < 1.5 with a constant initial ve- 
locity field v = [1, 1]. The initial gas pressure is P = 6 with 
7 = 5/3 and the magnetic field has a constant component 
perpendicular to the plane B z = l/\/4~7r. The divergence is 
introduced as a peak in the x— component of the field in the 
form 

B x = (r/r ) 8 - 2(r/r,,) 4 + 1 r = sj x 2 + y 2 (84) 

where ro is the radius of t he initial peak . The s etup used here 
differs from that used bv lDedner et al.| (I2002T) in that ro is a 
changeable parameter (using ro = l/VS gives their setup). 
The reason for this is that we find that the effectiveness 
of the divergence cleaning strongly depends on the wave- 
length of the divergence errors. Testing divergence clean- 
ing procedures bas ed on a single wavelength of error (as in 
iDedner et al"l l2002) can lead to misleading interpretations 
and an incorrect choice of parameters when applied to di- 
vergence errors which are generated in the course of real 
simulations. 

The contours of the initial B x field arrangement (1841 in 
the lDedner et all (|2002fl case of r = 1/V& are shown in the 
left hand side of Figure (and similarly in Figure |3Jl The 
particles are arranged on a cubic lattice for simplicity and in 
the evolution calculations the periodic boundary conditions 
are enforced using ghost particles. Since the density is uni- 
form throughout the simulation the results are insensitive 
to whether B or B/p is evolved and also to the instability 
correction method since the simulation is not unstable to 
negative stress. The artificial dissipation terms are turned 
off for this problem in order to isolate the effects of the di- 
vergence cleaning procedures. 

7.2.1 Projection methods 

We use this problem to test the projection methods by ap- 
plying a single projection step to the divergence error intro- 
duced in the initial conditions. To illustrate the divergence- 
free configuration adopted by the field we plot in Figure 
the contours of B x before and after the projection step. The 
initial configuration (left panel) is set according to 1841 with 
rp = l/y/E, corresponding to that used by IDedner et alJ 
(2002). The right panel shows the resulting field configura- 
tion after the projection step. We have used 2,500 (50 x 50) 
particles in this case, setup on a cubic lattice with a smooth- 
ing length set using r/ = 1.2 in 1131 . giving h — 0.048. This 
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Figure 1. Divergence cleaning using the approximate (scalar) 
projection method described in Wj.L'l The plot shows 30 contours 
of B x in the V-B advection problem before (left) and after (right) 
a single projection step at t = 0. The projected magnetic field 
adopts an essentially divergence- free configuration in a single step. 
Note that the wavelength of the initial divergence error in this 
case is substantially larger than the smoothing length. 



scalar projection 

vector projection 



Figure 2. Wavelength dependence of the scalar and vector ap- 
proximate projection methods described in i!6.2l The y— axis 
shows the relative change in the maximum divergence error (ie. 
max(V-B new /V-Bo) over a single projection step taken at t = 0. 
The x— axis gives the radius of the initial peak in the x component 
of the field ro in units of the smoothing length h. In both cases 
the projection step becomes less effective as the wavelength of the 
divergence error approaches the smoothing length (ie. ro/h — > 1), 
however the vector projection outperforms the scalar projection 
by a factor of ~ 100 at all wavelengths. 



means that the wavelength of the initial divergence error 
is substantially (~ 7x) larger than the smoothing length. 
However, we expect that divergence errors generated in the 
course of real simulations will tend to have wavelengths 
closer to ~ h. For this r eason we have extended the test 
problem of lDedner et al.l J2002t) to a variety of wavelengths 
by varying the parameter ro. 

The relative divergence cleaning given from a single pro- 
jection step using both the scalar ( H6.2.H and vector ( H6.2.2H 
projection methods are plotted against ro/h in Figure|2] The 
particles have been setup as previously, in this case varying 
ro whilst keeping h fixed. We have also performed simu- 
lations where ro is fixed and h is varied by changing the 
number of particles. The results in both cases are virtually 
identical. The results are also insensitive to the absolute size 



of the divergence error, since we plot only the relative change 
in V ■ B. 

For simplicity we have assumed that the boundaries are 
open when calculating the sum in the projection step, rather 
than explicitly treating the periodic boundaries. This is a 
reasonable assumption whilst the source terms for the Pois- 
son equation (ie. V • B or V x B) are non-zero in only a 
finite region of the simulation volume. However, to ensure 
that the maximum V ■ B value is not an artefact of our 
non-treatment of the periodic boundary conditions, the di- 
vergence near the boundaries of the domain (on particles 
within 2h of the boundary) has been set to zero when cal- 
culating the maximum used in Figure [5] These calculations 
shown in Figure |2] have also been performed with the B x 
peak placed in the centre of the domain rather than at the 
origin in order to move the source further away from the 
boundaries. 

It can be seen from Figure [5] that the effectiveness of 
the divergence cleaning given by the scalar projection is re- 
duced as the wavelength of the divergence error approaches 
the smoothing length (for the scalar projection the reduc- 
tion in error approaches a mere factor of ~ 2 as ro/h — * 1). 
Whilst the vector projection (dashed line) shows a similar 
trend with wavelength, it is clear that this method is a vast 
improvement over the scalar projection, reducing the max- 
imum divergence error by a factor of ~ 100 compared to 
the scalar case. This is due to the fact that the vector pro- 
jection gives an expression for B which analytically has a 
zero divergence (refer to the discussion in H6.2.21 . The only 
non-zero divergence resulting in this case is due to the SPH 
operator used to compute V • B after the projection step. 
This is therefore a substantial improvement over the scalar 
projection, in which V ■ B is not guaranteed to be exactly 
zero in any approximation. We have also experimented with 
the scalar projection using SPH operators for V • B other 
than all of which we find show similar results. 



7.2.2 Source terms and Hyperbolic/parabolic cleaning 

The source term approach ( H6.1H and the hyper- 
bolic/parabolic cleaning ( Hfcj.31 are tested by evolving the 
initial field configuration given by 1841 forward in time. The 
results of this test are shown in Figure^] The plots show the 
divergence of the magnetic field as it evolves in each case. 
The results using the consistent formulation of V • B terms 
discussed in paper I and in Htj.ll are shown in the top row. 
In this case the divergence error is passively advected by 
the flow and both the field and the divergence error remain 
unchanged (relative to the flow) at t = 1, demonstrating 
that the formalism is indeed consistent in the presence of 
magnetic monopoles and conserves the integral 1551 . 

In order to compare these results with a conservative 
[in the sense of conserving E5H formulation of the MHD 
equations, we have performed a simulation using an SPH 
induction equation of the form 



d 

Tit 



E 



nib 



Pa 



B 3 



dWai 

dxi 



(85) 



which is an SPH form of the conservative (in a volume sense) 
induction equation 
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Figure 3. Results of the V ■ B advection problem. An initially non-zero divergence is setup as a peak in the a:— component of the 
magnetic field (leftmost figures), with a velocity field v(x, y) = [1, 1] and periodic boundaries. The plots show renderings of V ■ B in the 
range — 1 < V ■ B < 1 (from black to white) at various times throughout the simulation for various divergence cleaning procedures. The 
consistent treatment of V • B terms (top row) is clearly seen to advect the divergence without change, which is an improvement over a 
"conservative" formulation of the MHD equations in which the divergence is smeared throughout the simulation volume (second row). 
With the use of hyperbolic cleaning in addition to the consistent V ■ B terms, the divergence error is spread rapidly in a wavelike manner 
(third row), whilst with a mixed hyperbolic/parabolic cleaning (fourth row) this error is also quickly diffused away. 



dt 



B 



B 



■ V v + v 



V B 



(86) 



The results using this formalism are shown in the second 
row of Figure|3] The peak in B x is distorted by the flow and 
the divergence error is smeared throughout the simulation. 

The third row in Figure |3 shows the results using the 
divergence correction discussed in H6.3l using only the hyper- 
bolic term in (1711 (ie. with a — 0) in conjunction with the 



usual monopole formulation of the induction equation (11911 . 
The divergence error is spread rapidly in a wavelike manner 
by the constraint equation. However, the magnitude does 
not decrease substantially in this case. 

Using the mixed hyperbolic/parabolic cleaning with a 
small amount of diffusion (using the parabolic term in 17111 . 
in this case with a — 0.1), this error is rapidly diffused away, 
resulting in a divergence-free field configuration (Figure [31 
bottom row). For comparison, the results of a single pro- 
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Figure 4. Time evolution of the maximum (left) and average 
(right) values of |V • B| over the particles. With a conservative 
formulation of the induction equation the divergence error in- 
creases with time (dashed line) whereas the errors are conserved 
using a formulation which is consistent in the presence of magnetic 
monopoles (solid line). With hyperbolic cleaning (dot-dashed) the 
maximum is quickly reduced although the average increases, how- 
ever with the parabolic term included the average error is also 
rapidly diffused away (dotted line). 



jection step at t = are shown in Figure Q showing the 
divergence-free configuration adopted by the field. 

The time evolution of various quantities throughout 
these simulations are shown in Figure |U The left panels 
show the evolution of the maximum (top) and average (bot- 
tom) of |V • B|. In conservative form (solid line) the maxi- 
mum divergence varies slightly and initially becomes larger 
than the initial value. The bottom panel shows that the av- 
erage value in this case steadily increases over time, due to 
the smearing effect of the divergence propagation 1501 . The 
consistent formulation of V • B terms (dashed line) main- 
tains a steady value of both the maximum and average, as 
observed in Figure^] With hyperbolic cleaning (dot-dashed) 
the maximum divergence error is quickly reduced (although 
increases at later times as the divergence waves cross the 
periodic domain and interact) whilst the average climbs as 
the divergence error is spread throughout the domain. Using 
the mixed hyperbolic/parabolic cleaning as described above 
(dotted line), both the maximum and average divergence is 
swiftly reduced. 

Finally it is important to examine the effect of varying 
the strength of the parabolic (diffusion) term in 1711 . The ef- 
fects of varying the diffusion parameter a for this particular 



Figure 5. Circularly polarized Alfven wave test. The left figure 
shows the particle setup in the lowest resolution simulation. On 
the right the vertical component of the magnetic field is plotted as 
a rendered image from the 32 X 64 particle run at t = 5, showing 
the propagation of the wave with respect to the domain and the 
particle setup. 



problem were explored in IPricd <l2004t) . It was later realised 
however that these results depended strongly on the value of 
the smoothing length (ie. the resolution of the calculations). 
Despite this the optimal parameter in other simulations was 
found to be independent of h. The reason for this is quite 
simple in hindsight and is due to our use of h as the length 
scale (wavelength of critical damping) in 175H . In this test 
problem the divergence error is setup in the initial condi- 
tions at a wavelength ro which is independent of the actual 
resolution used. Thus using a length scale h in 1751 . the op- 
timal cleaning is strongly dependent on resolution, since the 
optimal wavelength in this case corresponds closely to the 
wav elength of the i nitial divergence error. This also explains 
whv lDedner et all (|2002) found that their use of a fixed pa- 
rameter cy which has dimensions of length was found to 
give cleaning which is independent of resolution (but only 
for this specific problem!). Since in realistic calculations di- 
vergence errors are produced at wavelengths ~ h, in general 
the length scale used in 1751 should reflect this. We there- 
fore retain the length scale h but defer examination of the 
effect of varying the parameter a to the Orszag-Tang vortex 
problem ( H7.7H where divergence errors are generated in the 
course of the evolution. 



7.3 Circularly polarized Alfven wave 

This test is described bv lTotbl ( 2000) where it is used to test 
a variety of multidimensional MHD schemes in grid based 
codes. The test involves a circularly polarized Alfven wave 
propagating in a two dimensional domain. The advantage of 
using a circularly (as opposed to linearly) polarized wave is 
that it turns out to be an exact, non-linear solution to the 
MHD equations, which means that the solution after one 
period should exactly match the initial conditions, without 
the effects of nonlinear steepening (as observed, for example, 
in the magnetosonic wave tests described in paper II). This 
also means that the wave can be setup with a much larger 
amplitude than would be used for purely linear waves. 
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Figure 6. Results of the circularly polarized Alfven wave test at t = 5 (corresponding to 5 wave periods). The plots show the perpendicular 
component of the magnetic field vector B± = B y cos 6 — B x sin 9 for all of the particles, projected against a vector parallel to the direction 
of wave propagation rn = X cos + j/ sin (where 8 = 30° in this case). The SPMHD results are shown at five different resolutions which 
are, from bottom to top, 8 X 16, 16 X 32, 32 X 64, 64 X 128 and 128 X 256. Initial conditions are indicated by the solid line. The numerical 
results should match these initial conditions at the time shown. The left panel shows the results in the absence of dissipative terms and 
demonstrates that the SPMHD algorithm contains very little intrinsic numerical dissipation even at low resolutions, although there is a 
small phase error present even in the converged higher resolution runs. The right hand panel shows the results applying the dissipative 
terms required in the shock tube problems uniformly (ie. in the absence of switches). In this case the wave amplitude is damped by the 
artificial resistivity term and exhibits somewhat slow convergence. 



In lTothl |2000T) . the wave is setup to propagate at an an- 
gle 8 = 30° with respect to the a;— axis. In SPH the orienta- 
tion of the wave vector with respect to the co-ordinates is not 
particularly important because there is no spatial grid. How- 
ever, we have retained the rotated configuration as firstly it 
ensures that there are no spurious effects resulting from the 
initial arrangement of the particles and s econd l y ena bles a 
fair comparison with the results shown in iTothlfeood) . The 
particles are setup on a hexagonal close packed lattice (ie. 
such that particles are equispaced) in a rectangular domain 
< x < 1/ cos#;0 < y < 1/ sin#. This positioning of the 
boundaries means that periodic boundary conditions can be 
used, although some care is required to ensure the continu- 
ity of the lattice across the boundaries. This is achieved by 
stretching the lattice slightly in the y— direction to ensure 
that the boundaries lie at exactly half the spacing of the rows 
in the lattice. The particle setup at the lowest resolution is 
shown in the left hand side of Figure 

The wave is setup with a unit wavelength along the 
direction of propagation (ie. in this case along the line 
at an angle of 30° with respect to the x-axis). The ini- 
tial conditions are p = 1, P = 0.1, tin = 0, Bii = 1, 
v± — B± — 0.1 sin (27rrii ) and v z — B z = 0.1 cos (27rri| ) with 
7 = 5/3 (where ry = x cos 8 + y sin 8). The x— and y— com- 
ponents of the magnetic field are therefore given by B x — 
Bit cos 8 — B± sin 8 and B y = Bu sin 8 + B± cos 8 (and simi- 



larly for the velocity). Conversely, Bu = B y sin 8 + B x cos 8 
and B± — B y cos8 — B x sin8. Note that this setup means 
that V • B = holds as a combination of the dB x /dx and 
dB y /dy terms, rather than both components being zero in- 
dividually. The vertical component of the magnetic field af- 
ter 5 periods is plotted as a rendered image in the right hand 
side of Figure El showing the direction of wave propagation 
with respect to the domain and the particle setup. 

We have peformed this test at five different resolutions: 
8 x 16, 16 x 32, 32 x 64, 64 x 128 and 128 x 256 particles. 
In each case the number of particles in the y-direction is de- 
termined by the hexagonal lattice arrangement. The results 
are shown in Figure|H]after 5 wave periods (corresponding to 
t = 5). The plots show the perpendicular component of the 
magnetic field B± plotted against rii for all of the particles 
in the simulation, with the results from the bottom to top 
panels shown in order of increasing resolution. In each case 
the initial conditions are indicated by the solid line which is 
identical to the exact solution at the time shown. 

The left hand side of Figure |S| shows the results in the 
absence of dissipative terms (that is with the artificial viscos- 
ity, resistivity and thermal conductivity turned off). In this 
case the amplitude agrees very well with the exact solution 
even at the lowest resolutions. This demonstrates that SPH 
has a very low intrinsic numerical dissipation (compare for 
example with the damping of the wave at lower resolutions 
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in the plots shown in lT6thll200d) . However there is a small 
phase error which remains even in the highest resolution 
run. This is similar to the phase error ob served in the one 
dimensional sound wave tests presented in lPricd i2004t) and 
in the one dimensional magnetosonic waves tests in paper I. 
In these cases the phase error was found to be essentially 
removed by accounting for the variable smoothing length 
terms (paper II). The results shown in Figure |S| incorporate 
the variable smoothing length terms, however in this case 
the phase error is not completely removed (although is still 
an improvement over the results using simple averages of 
the smoothing lengths or kernel gradients) unless a number 
of neighbours used is also increased in addition to the total 
number of particles. 

The right hand side of Figure |H] shows the results of 
this test using the dissipative terms as required in the shock 
tube problems. In this case the wave is severely damped and 
convergence of the amplitude towards the exact solution is 
quite slow. The damping is largely caused by the uniform 
application of artificial resistivity (ie. using olb — 1 every- 
where) resulting in a somewhat large dissipation even in the 
absence of shocks. Substantially improved results could be 
obtained using the resistivity switch discussed in H5.2I how- 
ever for the shock tube problems it was found that use of 
such a switch could result in too little dissipation at rota- 
tional discontinuities in the absence of a shear viscosity term. 
Nonetheless the results shown in Figure||j]suggest that some 
kind of resistivity switch would be very valuable in SPMHD 
calculations. Note that the divergence error remains very 
small [(V ■ B) max ~ 10 -3 ] in all of the simulations shown 
even in the absence of any kind of divergence cleaning. 

7.4 2.5D shock tube 

The next two tests are simply two dimensional versions 
of the one di mensional sh ock tube tests described in pa- 
per I (see also lPricdl20o3) and demonstrate the effects of 
divergence errors in the shock capturing scheme. In two di- 
mensions we setup the particles on a cubic lattice in the 
a;— direction in the domain x = [—0.5 — v x (z,\t ma x, 0-5 — 
Vx(n)tmax\, where iWm and Vxim are the initial velocities 
assigned to the left and right states. This means that at 
the time t ma , x the particles are contained in the domain 
x = [—0.5,0.5]. The domain has a width of 4 particle spac- 
ings in the y— direction for computational efficiency. Bound- 
ary conditions are implemented by fixing the particle prop- 
erties in two buffer regions at the edges of the a;— domain, 
in which particles are evolved with a fixed velocity but copy 
their properties (p, P, B) from the nearest 'active' particle. 
Periodic boundary conditions are used in the y— direction, 
implemented using ghost particles. The exact position of 
the y— boundary is chosen to ensure periodicity of the lat- 
tice arrangement, ie. at half the spacing of the initial rows 
of particles in the y-direction. The initial shock is setup as 
a discontinuity in the fluid quantities at x = to which no 
smoothing is applied. 

The first shock test is the adiabatic shock tube problem 
involving seven different discontinuities given in paper I. 
Strictly this is a '2^' dimensional problem since the trans- 
verse velocity and magnetic field also have components in 
the z— direction. Conditions to the left of the discontinuity 
(the left state) are given by (p, P, v x , v y , v z , B y , B z ) = 



[1.08, 0.95, 1.2, 0.01, 0.5, 3.6/(4tt) 1/2 , 2/(4tt) 1/2 ] whilst 
to the right (the right state) the conditions are 
(p,P,v x ,v y ,v z ,B y ,B z ) = [l,l,0,0,0,4/(4 7 r) 1 / 2 ,2/(4^) 1 / 2 ] 
with B x = 2/(47r) 1 / 2 everywhere and 7 = 5/3. The problem 
has been studied by in one dimen sion by many authors (e.g. 
Rvu fc Joneslll995l:lBalsaralll998D and in two dimensions by 
T6thH2000h and lDedner et alJ J2003) . 

The problem was computed using 310x4 particles which 
corresponds to particles being uniformly spaced on a cubic 
lattice with separation 0.004 (with a slightly larger spacing 
for x > to give the density contrast), although results 
are similar using a hexagonal close packed lattice arrange- 
ment. Note that the above figure refers to the number of 
particles in the domain —0.5 < x < 0.5 at t max = 0.2 and 
that the resolution in this domain is correspondingly lower 
at earlier times due to the inflow boundary condition. The 
resoluti on wa s chos en to be comparable to the resolutions 
used in iTothl (|2000). The small density difference between 
the left and right states was setup by changing the lattice 
spacing slightly in the a;— direction. 

The solution using an initial smoothing length of h — 
1.2(m/p) 1//2 is shown in Figure|7]at t max = 0. 2 and may be 
compa red with the exact solution taken from lRvu fc Jonesl 
( 1995) (solid lin e) and with the one dimensional SPMHD 
results shown in lPricei i2004T) and in paper I. 

On the shock tests the most important physical aspects 
for a numerical algorithm are obtaining the correct physi- 
cal states behind the shock fronts, since this represents the 
manner in which the gas is changed by the passage of the 
shock. Thus whilst the shock profiles shown in Figure Q arc 
not as s harp as those shown at comparable resolution in, for 
example ITothl (2000), the intermediate states are obtained 
correctly apart from some small oscillations observed in the 
transverse velocity components near the contact discontinu- 
ity and the very narrow spikes in the magnetic field and 
transverse velocities which are damped at this resolution by 
our use of artificial resistivity. 

Whilst the damping due to artificial resistivity improves 
with resolution (and with the use of the resistivity switch - 
see Figure |5J, the oscillations in transverse velocity are of 
more concern. It should be noted first of all that these oscil- 
lations are quite small and do not appear to affect the dy- 
namics significantly (mainly because the jumps in the trans- 
verse velocity components are an order of magnitude less 
than the jump in v x ). However, the oscillations appear to 
result from a combination of three factors: the unsmoothed 
initial conditions, the fact that we do not explicitly apply 
any smoothing to the transverse velocity components and 
the effects of the small jumps in the a;— component of the 
magnetic field. 

To remove these oscillations two approaches can be 
taken: The first approach is to modify the artificial viscosity 
terms slightly in order to smooth the transverse velocity pro- 
files. The dissipative terms used in order to capture shocks 
were discussed at length in paper I and in this paper in 33 
In the one dimensional case the dissipation terms for MHD 
(comprising an artificial viscosity, artificial thermal conduc- 
tivity and artificial resistivity) were derived assuming that 
jumps would only occur in components of the magnetic field 
transverse to the line joining the particles that jumps in ve- 
locity would only occur parallel to this line. Neither of these 
assumptions strictly hold in the shock tube problem shown 
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Figure 7. Results of the 2.5D shock tube test using 310 X 4 particles and an initial smoothing length of h = 1.2(m/p) 1//2 . In two 
dimensions at this value of smoothing length small oscillations in the transverse velocity components appear primarily as a result of the 
non-zero magnetic divergence. In this plot the usual artificial viscosity and resistivity terms have been applied uniformly (ie. not using 
switches). A small amount of artificial thermal conductivity has been applied using the switch. 






Figure 8. The parallel component of the magnetic field in the 2.5D shock tube problem using the usual dissipative terms (left), using 
the total magnetic energy (centre) and using the total magnetic and kinetic energies (right). Using the total magnetic energy in the 
dissipative terms means that jumps in the parallel field components are smoothed in addition to the jumps in transverse field. Using the 
total kinetic energy smooths jumps in the transverse (as well as parallel) velocity components, however this explicitly adds an undesirable 
shear component to the artificial viscosity term. Details of these formalisms are given in JS] 
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Figure 9. Results of the 2.5D shock tube test using 310 X 4 particles and an initial smoothing length of h = 1.5(m/p) 1//2 and the total 
magnetic energy in the artificial resistivity term but using the usual artificial viscosity term (where in this case both have been applied 
using the dissipation switches). The results are a substantial improvement on those presented in Figure Q for a very modest increase in 
the number of neighbours. 



in Figure [7| since the transverse velocity components clearly 
jump and there is also a small jump in the parallel field 
component due to the divergence errors. 

A reformulation of the dissipative terms relaxing both 
of these assumptions was presented in H5.ll deriving the ar- 
tificial viscosity and artificial resistivity terms from jumps 
in the total kinetic and magnetic energies respectively in 
the total energy equation. The effects of using these for- 
mulations on the profile of the parallel component of the 
magnetic field are shown in Figure [S] From the centre panel 
we see that using the total magnetic energy formulation for 
the artificial resistivity has clear advantages in preventing 
oscillations in the parallel component of the field at shock 
fronts. Using the total kinetic energy version of the artificial 
viscosity (in order to smooth out jumps in the transverse 
velocity) effectively adds an explicit shear component to the 
viscosity term. 

In H5.ll it was noted that discontinuities in the trans- 
verse velocity components can only occur at corresponding 
jumps in the magnetic field and therefore that such disconti- 



nuities are already smoothed somewhat by the application of 
artificial resistivity there. For this reason the total kinetic en- 
ergy formalism is not strictly necessary provided that there 
is sufficient artificial resistivity present to smooth both the 
transverse field jumps and the transverse velocity jumps. 
However, applying even a small amount of such a viscosity 
to the two dimensional problem is i ndeed found to remove 
the observed oscillations JPricel2004) . It is clear though that 
the use of this term is highly undesirable since applying an 
explicit shear viscosity will substantially increase the spuri- 
ous transport of angular momentum caused by the artificial 
viscosity term. 

The second approach is to simply increase the number 
of neighbours slightly for each particle to give a more ac- 
curate interpolation. The results using an initial smoothing 
length of h = 1.5(m/ p) 1 ^ 2 are shown in Figure [5] using the 
total magnetic energy formulation of the artificial resistiv- 
ity but retaining the usual artificial viscosity formulation. In 
this case the jump in the parallel field component is much 
lower and the oscillations in the transverse velocity com- 
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ponents do not appear, although there is a small glitch at 
the contact discon tinuity simil ar to that observed in the one 
dimensional case (|Pricel l2004f) . The increase in smoothing 
length also means that the dissipative terms can be applied 
using the switches discussed in H5.2I resulting in a much 
lower dissipation rate away from the shocks than would be 
required for the h = 1.2(m/p) 1//2 case. 

Increasing the smoothing length from h = 1.2(m/p) 1 ^ 2 
to h = 1.5(m/p) 1//2 corresponds to an increase in the num- 
ber of neighbours from « 20 to w 28 on a uniform cubic 
lattice in two dimensions. This quite a small increase in 
computational expense for a substantial gain in accuracy 
(and stability). It therefore seems much more desirable to 
increase the smoothing length slightly for multidimensional 
problems rather than to explicit add a shear viscosity term. 

Finally, although this problem is not unstable to the 
clumping instability (and indeed no clumping is observed) 
we have also investigated the effects of various instability 
correction methods on the shock profile. In particular use 
of the antidumping term (paper I) was found to produce 
additional noise in the shock profile. Using either the Mor- 
ris formalism for the anisotropic force ( H4.21 or subtracting 
the constant component of the magnetic field ( H4.ll both 
give results very similar to those shown in Figures 17191 Ap- 
plying the hyperbolic/parabolic divergence cleaning to this 
problem gives a small reduction in the divergence error but 
otherwise has no significant effect on the shock profiles. 

7.5 Two dimensional shock tube 

The second shock tu be test is used by both lTothl l)200(t) and 
iDedner et alJ £2002) in two dimensions to compare the re- 
sults of various divergence cleaning schemes, although the 
one dimensio nal version of this test h a s been used by man y 
authors (e.g iDai fc Woodward! Il994t LRvu fc Jonesl 11995). 
The results of the one dim ensional test using the SPMHD 
algorithm are presented in IPricel l|2004h . Although this is 
a purely two dimensional test we present it after the 2.5D 
shock tube since it presents a much more challenging prob- 
lem with regards to the non-zero divergence of the magnetic 
field due to the stronger shocks. 

The particle setup is as described in the previous 
section, except that the initial left state is given by 
(p,P,v x ,v y ,B y ) = [1,20, 10,0, 5/(4tt) 1/2 ] and the right state 
is {p,P,v x ,v y ,B y ) = [1,1,-10,0,5/(4tt) 1 / 2 ] with B x = 
5.0/(47r) 1 ^ 2 and 7 = 5/3. The boundaries are correspond- 
ingly adjusted in the a;— direction to allow the particles to fill 
the domain —0.5 < x < 0.5 at t m „ = 0.08. Particles are ar- 
ranged initially on a hexagonal lattice with particle spacing 
0.004, giving 660 particles in the x— direction and a total par- 
ticle number of 660 x 4 = 2640. As in the previous test, the 
results using an initial smoothing length of h = 1.2(m/p) 1//2 
exhibit significant oscillations in the transverse velocity (v y ). 
In this case the oscillations are substantially worse because 
the jump in the parallel field component is much larger. 
Hence we have performed this test using h = 1.5(m/p) 1//2 . 
However, even in this case the oscillations remain present 
and so we have also added the shear viscosity term, using 
1371 with a = 1 everywhere (that is, not using the viscosity 
switch). The results using these settings are shown in Fig- 
ure | 10l an d may be compared with the exact solution taken 
from lDai fc Woodward! ill 9941) (solid line) and with the one 



dimensional results given in IPricd i2004l) . All particles are 
shown projected along the x— direction. Even in this case 
some oscillations are visible in the v y profile, corresponding 
exactly with a spike in V ■ B. In the h = 1.2(m/p) 1//2 case 
this spike is much larger [(V • ~B) max ~ 40], causing signifi- 
cantly more disruption to the velocity profile. Thus despite 
the various tweaks we have attempted for this test, the os- 
cillations appear to be primarily caused by the divergence 
errors generated at the shocks. More importantly a slight 
offset in the intermediate states in the v y , B y , p and P pro- 
files is present. Investigation of this effect suggests that this 
is caused by a combination of the divergence error and the 
shear viscosity term. Without the shear viscosity term the 
intermediate states are obtained correctly, although with 
substantially more oscillations in the v y profile. 

The effects of increasing the number of neighbours and 
changing the strength of the dissipation terms may be sum- 
marised as follows: Increasing the number of neigbours re- 
duces the jumps in the parallel field component (for exam- 
ple with h — 1.2(m/p) 1 ^ 2 the jump is given by AB X = 
[B x ( max ) - B x(min) ]/B x0 ~ 18% whilst for h = 1.5(ra/p) 1/2 
we have AB X « 1% and for h = 2.4(m/p) 1//2 this reduces 
further still to AB X « 0.15%). On the other hand, adding 
dissipation at the jumps in parallel field means that although 
such jumps may be present, the discontinuities (causing 
strong divergence errors) are smoothed. The effect of adding 
the shear viscosity term is to increase the dissipation at these 
discontinuities, thus reducing to some extent the associated 
spike i n the magnet ic divergence. 

In iTothl J2000J) the results of this test were pres ented 
using the source term approach of lPowell et al.l lll999T) (dis- 
cussed in H6.ll . showing similar jumps in the parallel mag- 
netic field component which were unchanged even in the 
converged numerical results. The fact that the jumps in par- 
allel field reduce with an increasing number of neighbours 
indicates that the SPMHD algorithm converges to the exact 
solution in the limit of h — » 00 an d N — > 00 where N is 
the number of particles. ITothl i2000l) attributes the errors in 
the parallel field components in the Powell method to the 
non- conservative source terms in the induction equation. We 
have also performed this simulation using the 'conservative' 
induction equation 1851 . however we find that the jumps 
in B x are not changed significantly by including the vV • B 
term (although contain substantially more numerical noise) . 

The shock tube tests presented above have been com- 
puted without using any form of divergence cleaning (other 
than the consistent formulation of the MHD equations in 
the presence of magnetic monopoles discussed in H6.ll . Thus 
a way of eliminating both the jumps in parallel field and 
the resulting oscillations in the transverse velocity compo- 
nents is to clean up the divergence error. Using the hyper- 
bolic/parabolic cleaning discussed in H6.3l is not particularly 
effective for this problem, since the divergence errors are 
propagated away from their source at the fastest wave speed 
which is similar to the rate at which they are created by the 
shocks. Thus the diffusion introduced by the parabolic term 
does not have time to eliminate the divergence error be- 
fore oscillations in the velocity components are produced. 
This is illustrated in Figure 1111 which shows the parallel 
field component and the divergence error after using this 
type of cleaning with a = 0.4 in the parabolic term (c.f. 
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Figure 10. Results of the two dimensional shock tube test at t = 0.08 using h = 1.5(m/p) 1 / 2 and the shear viscosity term. The results 
may be compared with the exact solution given by the solid line. In this stronger shock tube problem the jumps in the parallel field can 
cause significant oscillations in the transverse velocity components due to the non-zero divergence terms. 





Figure 11. Parallel magnetic field (left) and the divergence er- 
ror (right) in the two dimensional shock tube test at t = 0.08 
computed as in Figure [Tul but using the hyperbolic/parabolic di- 
vergence cleaning f i|6.3l . The exact solution is given by the solid 
line. The hyperbolic divergence cleaning does not have a large 
effect on this problem since the divergence errors are propagated 
at the fastest wave speed which is similar to the rate at which 
they are generated in the shocks. 



The diver gence errors are reduced by a factor of ~ 2 
compared to the results shown in Figure tTUl In order to elim- 
inate the divergence errors from problems such as this one 
where divergence errors are created rapidly it seems neces- 
sary to invoke some kind of sub-timestep cleaning (such as 
a projection method) . The implementation of such methods 
are complicated in this simple test problem by the use of 
periodic boundary conditions. Alternatively the number of 
neighbours can be increased further. To demonstrate this 
we present a simulation at double the usual neighbour num- 
ber, that is using h = 2A(m/p) 1 ^ 2 . The results are shown 
at t = 0.08 in Figure 1121 and show a reduction in the di- 
vergence error by a factor of ~ 10 compared to the results 
shown in Figure [TUl This suggests that using a larger neigh- 
bour number may be crucial in three-dimensional SPMHD 
simulations. 



7.6 Rotor 

The next test is taken from iTothl (|200(J) and consists of a 
spinning, dense disc embedded in an ambient background 
medium containing a uniform magnetic field. The material 
initially contained within the disc is flung into the surround- 
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Figure 12. Results of the two dimensional shock tube test at t = 0.08 computed as in Figure [TUI but using h = 2.4(m/p) 1//2 . The exact 
solution is given by the solid line. Increasing the neighbour number significantly decreases the divergence error. 



ing medium by the centrifugal forces, but is constrained into 
an oblate shape by the magnetic field. The computational 
domain is given by —0.5 < x,y < 0.5 with uniform ther- 
mal pressure P = 1 and an adiabatic index of 7 = 1.4. A 
constant, uniform magnetic field is setup in the a;— direction 
with strength B x = The dense disc is setup with 

p = 10 and a rotation velocity given by v x = 2(y — 0.5)/ro, 
v y — 2(x — 0.5)/ro for r < ro where in this case ro = 0.1. 
The ambient medium is at rest with p = 1.0. Note that this 
choice of i nitial conditi ons corresponds to the 'first rotor 
problem' in lTothl feOOOt) . 

The density contrast between the disc and the back- 
ground medium can be setup in SPH using either variable 
particle masses and therefore a fixed initial separation or 
equal mass particles and a variable particle distribution. We 
have experimented with both methods. In the variable par- 
ticle mass case the large density contrast results in some 
spurious effects from the higher mass particles 'mixing' into 
the low particle mass medium and we therefore prefer the 
equal mass particle approach. We achieve this setup by set- 
ting up the initial disc with a dense concentration of particles 
setup on a regular, hexagonal close-packed lattice trimmed 
to r < ro. The surrounding medium is then placed using 
a second close packed lattice with a correspondingly larger 



inter-particle separation with the region r < ro excluded. 
This setup means that we do not apply a tap er function t o 
the density, pressure or velocity profiles as in iTothl <200Cl) . 
However the density profile is naturally tapered by the it- 
erative calculation of the smoothing lengths and densities 
of the particles across the interface (|3J|. To ensure numeri- 
cal pressure equilibrium we setup the thermal energy of the 
particles using uo = Po/[(j — l)po] after the initial density 
has been calculated by direct summation, rather than using 
the analytic density step. Despite this there are some initial 
transients but these do not appear to affect the subsequent 
evolution substantially. 

The problem has been calculated using a background 
medium with 200 particles in the x— direction. The hexag- 
onal lattice arrangement means that this corresponds to 
200 x 230 particles in the surrounding medium, from which 
the central disc region is removed, leaving 44,332 particles in 
the background medium. The dense concentration of parti- 
cles in the disc contains a further 14,454 particles, resulting 
in a total of 58,786 particles. Artificial viscosity and resis- 
tivity have been applied using the switches, with artificial 
thermal conductivity turned off. No divergence cleaning has 
been applied. 

The results at this resolution using h = \.2{m/ p)z 
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Figure 13. The density, pressure, mach number and magnetic pressure at t = .15 in the M HD rotor problem using 58,786 particles 
(roughly 200 X 200). All plots show 30 contours spaced between the limits given in lTothl <2000]) . that is 0.483 < p < 12.95, 0.0202 < P < 
2.008, < \v\/c a < 1-09, < |B 2 < 2.642 in order to make a direct comparison. 



are plotted in Figure 1131 and may be directly com pared 
with the high resolution results shown in iTothl teOOd) . The 
density resolution in the SPMHD solution is slightly bet- 
ter than ev en the 400 x 400 grid based solution shown in 
ITothl (2000), giving a maximum density of p m ax = 15.54 at 
t — 0.15 as opposed to pmax = 12.95 in the grid solution, al- 
though the minimum density at this time IS pmin — 0.74 
in th e SPM HD solution as opposed to p m in = 0.483 in 
ITothl fcOOOl) . The SPMHD result using 400 particles in the 
a;— direction (giving a total of 235,574 particles) resolves a 
density range of p ma x = 17.76 and p m in = 0.58. The max- 
imum field strength is a little lower in the SPMHD calcu- 
lations, with (\B 2 )rn,ax = 2.3 (or 2.45 at the higher resolu- 
tion), as opposed to (^B 2 ) max = 2.64 in the 400 x 400 grid 
solution. This is due to our use of artificial resistivity for 
shock capturing. There are some small effects at low densi- 
ties in the SPMHD solution due to the particle distribution. 
These effects decrease both with particle number and also as 
the number of neighbours is increased. The divergence con- 



straint is maintained reasonably well in this problem - for 
example 95% of the particles have /t|V-B|/|B| < 0.01 in the 
200 x 200 particle simulation, which increases to 98% using 
400 x 400 particles and decreases to 87% using 100 x 100 
particles. 

7.7 Orszag-Tang vortex 

The final test is the compressible Orsz ag-Tang vortex 
proble m which was first investigated by lOrszag fc Tand 
(1979) in order to study incompressible MHD turbu- 
lence. The problem was later exte nded to the com pressible 
case b y iDahlburg fc Piconel dl989h and iPicone fc Dahlburel 
dl99ll) . More recently it has been widely used as a 
test probl e m for multidimensional MH D algor ithm s (e.g . 
Rvu et alJ Il995t iBalsaxal Il998t jDai fc Woodward! Il998i : 
Londrillo fc Del Zannall2000l : lT6thll2000l) . 

The setup consists of an initially uniform density, pe- 
riodic lxl box given an initial velocity perturbation 
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Figure 14. Results of the two dimensional Orszag-Tang vortex test, showing the density (top) and magnetic pressure (bottom) distri- 
butions at t = 0.5 for resolutions of 128 X 146 (left), 256 X 294(centre) and 512 X 590(right) particles. The particles are initially arranged 
on an isotropic hexagonal lattice with periodic boundary conditions. The initial velocity field is a large vortex v = [— sin (2ixy), sin (2itx)} 
whilst the magnetic field has a doubly periodic geo metry B = Bg\— sin (27ry) , sin (47rx)l . The SPMHD results at higher resolutions are 
in excellent agreement with those presented in (e.g.) lDai fc Woodward 11998T) and iTothl (2000). 



v = vo[— sin (2iry), sin (2ttx)] where vo — 1. The mag- 
netic field is given a doubly periodic geometry B = 
Bo[— sin (27ry), sin (47rx)] where Bo = 1/s/4tt. The flow has 
an initial average Mach number of unity, a ratio of mag- 
netic to thermal pressure of 10/3 and we use 7 = 5/3. 
The initial gas state is therefore P — 5/3-Bq = 5/(12-7r) 
and p = 'yP/vo = 25/(367t). Note that the choice of length 
and time scales differs slightly between various implementa- 
tions in the litera ture. The setup used above follows that of 
iRvu et alJ l)l995h and lLondrilio fc Del Zannal (|2000h . 

The particles are arranged initially on a uniform hexag- 
onal close packed lattice. This distribution means that the 
particle are isotropically arranged and is the distribution 
towards which other arrangements naturally settle. How- 
ever, results are similar using a cubic lattice arrangement. 
The simulation is performed at three different resolutions: 
128 x 146, 256 x 294 and 512 X 590 particles (where the 
number of particles in the y— direction is determined by the 
isotropic lattice arrangement). The periodic boundary con- 
ditions are implemented using ghost parti cles. These resolu- 
tions are similar to the resolutions used in lDai fc Woodward! 



( 1998) (although in SPH the resolution is concentrated pref- 
erentially towards regions of high density). The dissipation 
terms are applied using the artificial viscosity and resistiv- 
ity switches but leaving the artificial thermal conductivity 
turned off in order to increase the density resolution. The 
wall heating effects which t he art i ficial thermal conductiv- 
ity prevents are discussed in lPricel l)2004F) and are in general 
quite minor. No shear viscosity term has been used. Simula- 
tions of this problem which have been run with or without 
the variable smoothing length terms, using the Morris for- 
malism for the magnetic force ( H4.21 . evolving either B or 
B/p and either the thermal or total energy show essentially 
no difference in the numerical results. 

The results of the density evolution are shown in Figure 
1141 at t = 0.5. At this time four main shock fronts are visi- 
ble which have interacted in the central regions after having 
crossed the periodic domain. The SPMHD results, particu- 
larly are in good agreement w ith those presen ted in (e.g.) 
iDai fc Woodward! (Il994l.ll998l) and lTothl fcOOfJ) . In the low- 
est resolution run, the central regions appear to be slightly 
better resolved than in the 128 x 128 fixed-grid simulation of 
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Figure 15. Effect of the hyperbolic/parabolic cleaning on the 
evolution of the average magnetic divergence in the two dimen- 
sional Orszag-Tang vortex problem, varying the parameter a. 
With cr too low the cleaning can cause increases in V • B (dashed 
line) over simulations with no divergence cleaning (solid line). 
The optimal cleaning is obtained with a ~ 0.4 — 0.8 (dot-dashed, 
dotted lines). However, the reduction in the divergence obtained 
using the hyperbolic cleaning is fairly small. The single biggest 
factor determining the magnitude of the divergence error is the 
number of neighbours. The results shown are for a smoothing 
length of h = 1.2(m/p) 1 / 2 , although the errors decrease as the 
number of neighbours is increased. 



iDai fc Woodward! jl99ct) , although the lower density regions 
are correspondingly less well resolved. At this resolution the 
SPMHD solution shows some small residual effects due to 
the distortion of the initial regular particle arrangement, 
noticable as small ripples behind the shock fronts in Figure 
1141 and a slightly mottled appearance in the low density re- 
gions. This is particularly evident in Figure [T4l since we have 
used a smoothing length of h = 1.2(m/p) 1//2 . In the lowest 
resolution run the density maxima visible in the higher reso- 
lution runs at the top and bottom of the domain are largely 
washed out. This is a result of the artificial resisitivity term 
used for shock capturing which dissipates energy in these 
regions due to the strong current gradient. 

The evolution of the average of the magnetic diver- 
gence is shown in Figure 1151 for several runs using the hy- 
perbolic divergence cleaning. The results using the hyper- 
bolic/parabolic cleaning with a = 0.2 (dashed line) can in 
fact increase the divergence error over the results with no 
divergence cleaning (solid line). This is because in the ab- 
sence of sufficient diffusion the hyperbolic term can spread 
the divergence errors such that the resultant 'divergence 
waves' can constructively interfere with each other, lead- 
ing to increased errors. The optimal cleaning is obtained 
with a ~ 0.4 — 0.8 (dot-dashed, dotted lines), although the 
reduction in the divergence error given by the hyperbolic 
cleaning is comparatively small. In fact, as in the previous 
tests, the single biggest factor which determines the magni- 
tude of the divergence error is the number of neighbouring 
particles. For example in a simulation using h = 1.5(m/p) 1//2 
the divergence errors are approximately half those shown in 
Figure ^| 



8 DISCUSSION 

In this paper multidimensional aspects of the SPMHD algo- 
rithm have been discussed. In particular several methods for 
maintaining the divergence-free constraint in an SPH con- 
te xt have been presen ted. Firstly the source term approach 
of lPowell et al.l l)l999f) was outlined and contrasted with the 
consistent formulation of the MHD (and SPMHD) equations 
derived in paper II. The major difference between the two 
approaches is that our approach ret ains the conse rvation of 
momentum and energy whereas the IPowell et al.l approach 
does not. The conservation properties of the induction equa- 
tion were also discussed, in which it was highlighted that 
using a 'non-conservative' induction equation means that 
the surface integral of the magnetic flux is conserved, rather 
than the volume integral. The effect of using the consistent 
formulation of the MHD equations in the presence of mag- 
netic monopoles (which conserves the surface integral of the 
flux) is that divergence errors are advected without change 
by the flow (illustrated in Figure 

Projection methods for maintaining a divergence free 
field were discussed in an SPH context in H6.2I In partic- 
ular it was noted that using the Green's function solution 
to the Poisson equation (as is often used for self-gravity in 
SPH) provides only an approximate projection. The results 
using this type of projection on a problem where an ini- 
tial magnetic divergence was introduced into the simulation 
were very good ( H7.2H . but were found to degrade as the 
wavelength of the divergence error approached the resolution 
length. A projection method based on Biot-Savart's law was 
also discussed and found to give excellent results even for 
wavelengths approaching the smoothing length. The imple- 
mentation of either of these projection schemes for the test 
problems considered in this paper was complicated by the 
periodic boundary conditions used, leaving a need for fur- 
ther testing of these methods on three dimensional problems. 
In particular the Biot-Savart projection method suggests a 
promising divergence cleaning method in three dimensions. 
This would however require implementation in a tree-code 
which is beyond the scope of this paper. 

An alternative approach to diverg ence cleaning sug- 
gested recently by iDedner et alJ l)2002f) was discussed in 
H6.3I The method involves adding an additional constraint 
equation which is coupled to the induction equation for the 
magnetic field. Chosen appropriately, the effect of this equa- 
tion is to cause the divergence errors to be propagated in a 
wave-like manner away from their source (Figure|2Jl. Adding 
a small diffusive term means that the divergence errors are 
also rapidly reduced to zero. This method is extremely sim- 
ple to implement and is computationally very inexpensive. 
The disadvantage is that the error propagation is limited by 
the timestep condition and, although much faster than using 
diffusion alone to reduce the divergence, for some problems 
(for example the shock tube tests given in HT.4I and 
the cleaning is still not fast enough. However, this method 
is some improvement over not using any form of divergence 
cleaning at a negligible additional computational cost. 

The SPMHD algorithm was also tested against a variety 
of multidimensional test problems. A non-linear circularly 
polarized Alfven wave was studied in H7.MI This test showed 
that SPMHD has a very low intrinsic numerical dissipation 
compared to grid based codes, although this property is de- 
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stroyed by the addition of explicitly dissipative terms for 
shock-capturing which can cause quite slow convergence on 
problems where the physical dissipation timescale is of crit- 
ical importance. 

Two of the shock tube problems ex amined prev iously 
in one dimensional simulations (paper I IPricel 12004') were 
examined in two dimensions in H7.4l and l7.5l For these prob- 
lems jumps in the component of the magnetic field parallel 
to the shock front (causing divergence errors) were found 
to result in oscillations in the transverse velocity profiles. 
The jumps in the parallel field component were found to 
decrease as the number of neighbours for each particle was 
increased. The corresponding divergence errors produced by 
these jumps could be reduced by using a form of the dissipa- 
tive terms derived in 85. II using the total jump in magnetic 
and kinetic energies. Modifying the artificial viscosity term 
in this manner results in the addition of an explicit shear 
viscosity component. It is therefore somewhat undesirable 
to do so since this can result in excess spurious angular mo- 
mentum transport elsewhere. A better approach would be to 
use divergence cleaning to prevent these errors from occur- 
ing. However, the hyperbolic cleaning was not found to be 
particularly effective for this problem because of the restric- 
tion to the fastest wave speed and implementation of the 
projection method is complicated by the periodic boundary 
conditions. These difficulties are not, however, insurmount- 
able. The single biggest factor in determining the magnitude 
of the divergence errors in the shock tube tests was found to 
be the size of the smoothing region (ie. the number of con- 
tributing neighbours). It therefore seems advantageous to 
use a slightly larger number of neighbours for MHD prob- 
lems (typically h > 1.5(m/p) 1 ^" where v is the number of 
spatial dimensions) than might otherwise be used for hydro- 
dynamics. 

An MHD rotor problem was examined in 87.6 | with re- 
sults comparable to those shown in iTothl feOOOh . Finally 
the algorithm was tested on the Orszag-Tang vortex prob- 
lem (" H7.7I which has been widely used as a benchmark for 
MHD codes. The SPMHD results were in good agreement 
with those presented elsewhere. This test again highlighted 
the need for a slightly larger number of neighbours, in this 
case to remove spurious effects related to the initial lat- 
tice arrangement and to reduce the magnitude of the di- 
vergence errors. The hyperbolic/parabolic divergence clean- 
ing was found to produce only a small reduction in the di- 
vergence errors, again highlighting the need for some form 
of sub-timestep cleaning (for example using the projection 
method) . 

An issue which has not been discussed in this paper, 
but which needs to be addressed elsewhere, is the tendency 
of SPH particles merge together at short separations due to 
the fact that the force tends to zero near the origin of the 
cubic spline kernel. In particular this problem can become 
more acute as the number of neighbours is increased (as 
is required in order to maintain the divergence constraint 
in MHD). This instability is well known but is not neces- 
sarily noticeable in SPH simulations, particularly in 3 di- 
mensio ns, as it simply leads to a low er effective resolution. 
Whilst iThomas fc Couchmanl il992T) propose a simple solu- 
tion whereby the kernel gradient in the cubic spline is modi- 
fied slightly whilst retaining the usual kernel for the density 
evaluation, it is not clear what effect this has on the evolu- 



tion, particularly when using the variable smoothing length 
formalism which we have described here and in paper II. 
We therefore feel that this problem in particular warrants 
further attention. 

Finally it is worth commenting on the ability of the 
algorithm as it stands to treat 'real' astrophysical MHD 
problems. The crucial issue here is the degree to which the 
divergence constraint can be maintained. Of the methods 
examined in this paper the most promising is the projection 
method using the Biot-Savart law since it is the only method 
which guarantees a zero divergence. Efficient implementa- 
tion of this method in three dimensions requires use of a tree 
code (or similar) to solve the resulting Poisson-type equa- 
tion similar to that used to compute the gravitational force. 
Note that the treecode implementation differs slightly from 
the usual gravity tree since the source term of the Poisson 
equation in this case is a vector quantity. Periodic bound- 
ary conditions add a further complication although again 
methods used for gravity can be easily adapted. Secondly 
the issue regarding particle merging discussed above needs 
to be addressed to be able to usefully increase the neigh- 
bour number. Thus, whilst many improvements could still 
be made to the algorithm, the results presented in this paper 
suggest that the method is ripe for application to problems 
of current theoretical interest, such as that of star formation. 
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